Математичне моделювання кінетики росту кристала гідроксиапатиту
Розроблено математичну модель кінетики росту кристалів гідроксиапатиту, яка дає змогу проаналізувати вплив зміни кількості речовини вихідних реагентів і температури синтезу на розміри частинок. За допомогою квантовохімічних розрахунків отримані масиви залежностей термодинамічних параметрів в функції...
Saved in:
| Published in: | Поверхность |
|---|---|
| Date: | 2010 |
| Main Authors: | , , |
| Format: | Article |
| Language: | Ukrainian |
| Published: |
Інститут хімії поверхні ім. О.О. Чуйка НАН України
2010
|
| Subjects: | |
| Online Access: | https://nasplib.isofts.kiev.ua/handle/123456789/39322 |
| Tags: |
Add Tag
No Tags, Be the first to tag this record!
|
| Journal Title: | Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| Cite this: | Математичне моделювання кінетики росту кристала гідроксиапатиту / А.П. Головань, І.В. Головань, Є.М. Дем'яненко // Поверхность. — 2010. — Вип. 2(17). — С. 51-62. — Бібліогр.: 19 назв. — укр. |
Institution
Digital Library of Periodicals of National Academy of Sciences of Ukraine| id |
nasplib_isofts_kiev_ua-123456789-39322 |
|---|---|
| record_format |
dspace |
| spelling |
Головань, А.П. Головань, І.В. Дем'яненко, Є.М. 2012-12-15T10:57:28Z 2012-12-15T10:57:28Z 2010 Математичне моделювання кінетики росту кристала гідроксиапатиту / А.П. Головань, І.В. Головань, Є.М. Дем'яненко // Поверхность. — 2010. — Вип. 2(17). — С. 51-62. — Бібліогр.: 19 назв. — укр. XXXX-0106 https://nasplib.isofts.kiev.ua/handle/123456789/39322 544.4:544.18 Розроблено математичну модель кінетики росту кристалів гідроксиапатиту, яка дає змогу проаналізувати вплив зміни кількості речовини вихідних реагентів і температури синтезу на розміри частинок. За допомогою квантовохімічних розрахунків отримані масиви залежностей термодинамічних параметрів в функції температури та виведені їх поліноми. Проведений чисельний експеримент та отримані характеристики динаміки зміни кількості частинок за розмірами для різних температур синтезу . Разработана математическая модель кинетики роста кристаллов гидроксиапатита, которая дает возможность проанализировать влияние изменения количества вещества исходных реагентов и температуры синтеза на размеры частиц. С помощью квантовохимических расчетов получены массивы зависимостей термодинамических параметров от температуры и выведены их полиномы. Проведен численный эксперимент и получены характеристики динамики изменения количества частиц по размерам для различных температур синтеза. A mathematical model of hydroxyapatite crystal growth kinetics which enables to analyse influence the change of initial reagents substance quantity and synthesis temperatures on particle sizes is developed. Arrays of dependences of thermodynamical parameters on the temperature are obtained with quantum-chemical calculations and their polynoms are deduced. Numerical experiment is carried out and characteristics of dynamics changes of quantity particles in the sizes for various synthesis temperatures are obtained. uk Інститут хімії поверхні ім. О.О. Чуйка НАН України Поверхность Теория химического строения и реакционной способности поверхности. Моделирование процессов на поверхности Математичне моделювання кінетики росту кристала гідроксиапатиту Математическое моделирование кинетики роста кристалла гидроксиапатита Mathematical modelling of the hydroxyapatite crystal growth kinetics Article published earlier |
| institution |
Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| collection |
DSpace DC |
| title |
Математичне моделювання кінетики росту кристала гідроксиапатиту |
| spellingShingle |
Математичне моделювання кінетики росту кристала гідроксиапатиту Головань, А.П. Головань, І.В. Дем'яненко, Є.М. Теория химического строения и реакционной способности поверхности. Моделирование процессов на поверхности |
| title_short |
Математичне моделювання кінетики росту кристала гідроксиапатиту |
| title_full |
Математичне моделювання кінетики росту кристала гідроксиапатиту |
| title_fullStr |
Математичне моделювання кінетики росту кристала гідроксиапатиту |
| title_full_unstemmed |
Математичне моделювання кінетики росту кристала гідроксиапатиту |
| title_sort |
математичне моделювання кінетики росту кристала гідроксиапатиту |
| author |
Головань, А.П. Головань, І.В. Дем'яненко, Є.М. |
| author_facet |
Головань, А.П. Головань, І.В. Дем'яненко, Є.М. |
| topic |
Теория химического строения и реакционной способности поверхности. Моделирование процессов на поверхности |
| topic_facet |
Теория химического строения и реакционной способности поверхности. Моделирование процессов на поверхности |
| publishDate |
2010 |
| language |
Ukrainian |
| container_title |
Поверхность |
| publisher |
Інститут хімії поверхні ім. О.О. Чуйка НАН України |
| format |
Article |
| title_alt |
Математическое моделирование кинетики роста кристалла гидроксиапатита Mathematical modelling of the hydroxyapatite crystal growth kinetics |
| description |
Розроблено математичну модель кінетики росту кристалів гідроксиапатиту, яка дає змогу проаналізувати вплив зміни кількості речовини вихідних реагентів і температури синтезу на розміри частинок. За допомогою квантовохімічних розрахунків отримані масиви залежностей термодинамічних параметрів в функції температури та виведені їх поліноми. Проведений чисельний експеримент та отримані характеристики динаміки зміни кількості частинок за розмірами для різних температур синтезу .
Разработана математическая модель кинетики роста кристаллов гидроксиапатита, которая дает возможность проанализировать влияние изменения количества вещества исходных реагентов и температуры синтеза на размеры частиц. С помощью квантовохимических расчетов получены массивы зависимостей термодинамических параметров от температуры и выведены их полиномы. Проведен численный эксперимент и получены характеристики динамики изменения количества частиц по размерам для различных температур синтеза.
A mathematical model of hydroxyapatite crystal growth kinetics which enables to analyse influence the change of initial reagents substance quantity and synthesis temperatures on particle sizes is developed. Arrays of dependences of thermodynamical parameters on the temperature are obtained with quantum-chemical calculations and their polynoms are deduced. Numerical experiment is carried out and characteristics of dynamics changes of quantity particles in the sizes for various synthesis temperatures are obtained.
|
| issn |
XXXX-0106 |
| url |
https://nasplib.isofts.kiev.ua/handle/123456789/39322 |
| citation_txt |
Математичне моделювання кінетики росту кристала гідроксиапатиту / А.П. Головань, І.В. Головань, Є.М. Дем'яненко // Поверхность. — 2010. — Вип. 2(17). — С. 51-62. — Бібліогр.: 19 назв. — укр. |
| work_keys_str_mv |
AT golovanʹap matematičnemodelûvannâkínetikirostukristalagídroksiapatitu AT golovanʹív matematičnemodelûvannâkínetikirostukristalagídroksiapatitu AT demânenkoêm matematičnemodelûvannâkínetikirostukristalagídroksiapatitu AT golovanʹap matematičeskoemodelirovaniekinetikirostakristallagidroksiapatita AT golovanʹív matematičeskoemodelirovaniekinetikirostakristallagidroksiapatita AT demânenkoêm matematičeskoemodelirovaniekinetikirostakristallagidroksiapatita AT golovanʹap mathematicalmodellingofthehydroxyapatitecrystalgrowthkinetics AT golovanʹív mathematicalmodellingofthehydroxyapatitecrystalgrowthkinetics AT demânenkoêm mathematicalmodellingofthehydroxyapatitecrystalgrowthkinetics |
| first_indexed |
2025-11-26T17:20:51Z |
| last_indexed |
2025-11-26T17:20:51Z |
| _version_ |
1850764202137878528 |
| fulltext |
Поверхность. 2010. Вып. 2(17). С. 51–62 51
УДК 544.4:544.18
МАТЕМАТИЧНЕ МОДЕЛЮВАННЯ КІНЕТИКИ РОСТУ
КРИСТАЛА ГІДРОКСИАПАТИТУ
А.П. Головань1, І.В. Головань2, Є.М. Дем'яненко2
1
Інститут хімії поверхні ім. О.О. Чуйка Національної академії наук України,
вул. Генерала Наумова, 17, Київ, 03164, Україна
2
Інститут електродинаміки, Національної академії наук України,
пр. Перемоги 56, Київ, 03680,Україна
Розроблено математичну модель кінетики росту кристалів гідроксиапатиту, яка дає
змогу проаналізувати вплив зміни кількості речовини вихідних реагентів і температури синтезу
на розміри частинок. За допомогою квантовохімічних розрахунків отримані масиви
залежностей термодинамічних параметрів в функції температури та виведені їх поліноми.
Проведений чисельний експеримент та отримані характеристики динаміки зміни кількості
частинок за розмірами для різних температур синтезу .
Вступ
Моделювання кристалічних структур і дослідження їх властивостей за
допомогою методів мінімізації енергії міжатомної взаємодії стає доступним завдяки
розробці програмного забезпечення і використання потужної обчислювальної техніки.
Спочатку задача пошуку стійкої структури за допомогою мінімізації її повної
енергії була вирішена для молекулярних органічних кристалів. Успіх приніс метод так
званих атом-атомних потенціалів [1]. За допомогою даного методу вдається розрахувати
геометрію упаковки молекул в структурі, енергію сублімації органічних кристалів і деякі
їх фізичні властивості. Той факт, що ця задача вирішена раніше за інші, пояснюється
наступними причинами: по-перше, взаємодія в таких кристалах описуються силами Ван-
дер-Ваальса з великим радіусом дії (між сусідніми молекулами) і, по-друге, вони
обмежені невеликим набором пар легких атомів, таких як водень, вуглець, азот, кисень і
деякі інші, з потенціалами взаємодії, що переносяться від кристалу одного складу до
кристалу іншого складу. Залишається тільки варіювати розташування однієї молекули
відносно іншої до тих пір, поки не буде знайдена стабільна конфігурація, що
відповідатиме мінімуму енергії. Більшість таких розв’язків може бути знайдено
порівняно простими обчислювальними засобами.
Значно складніша задача пошуку оптимальної структури у випадку неорганічних
кристалів і мінералів. В таких кристалах ван-дер-ваальсові сили виявляються лише
відносно малою добавкою в загальну взаємодію, яка забезпечується більш міцним іон-
іонним та ковалентним чи більш складним по характеру хімічним зв’язком. Крім того,
різноманітність контактів між сусідніми атомами є набагато більшою, ніж у випадку
органічних кристалів. Тому складно чи навіть просто неможливо створити
універсальний набір парних потенціалів, які були б придатні для широкого кола таких
об’єктів. В додаток до цього необхідно відмітити, що взаємодія між атомами в таких
кристалах не обмежується тільки найближчими сусідніми, але й розповсюджується на
багато більші відстані, що захоплюють весь кристал в цілому. З цих причин
універсальний підхід до всіх типів неорганічних кристалів і мінералів є поки що
неможливим.
Незважаючи на ряд труднощів, пов’язаних з дослідженням кристалічних систем,
зроблена спроба дослідити механізм росту кристалів гідроксиапатиту (ГАП) та
формування необхідних параметрів з метою подальшого їх використання в рівнянні
52
кінетики та оцінки достовірності таких результатів. Критерієм адекватності отриманих
розрахункових даних фізико-хімічним процесам, при виборі комп’ютерних і мате-
матичних методів дослідження, стане їх співставлення з експериментальними даними.
Об’єкти і методи дослідження
Для дослідження комп’ютерними методами кристала гідроксиапатиту вибрано
структуру більшості його природних зразків, що відповідає гексагональній будові.
Згідно з експериментальними даними [2] елементарна комірка гідроксиапатиту має
розміри а = 9,52 Å, c = 6,874 Å, та відноситься до просторової групи Р63/m, із
структурними параметрами, наведеними в табл. 1. Така структура кристала (рис. 1) була
реалізована в програмі для візуалізації та дослідження кристалічних структур –
Diamond 3.0.
Таблиця 1. Структурні параметри гексагональної структури гідроксиапатиту
Атом x y z
Ca1 1/3 2/3 0,00300
Ca2 0,25100 1,00000 0,25
P1 0,40000 0,36900 0,25200
O1 0,33300 0,48700 0,25300
O2 0,58800 0,46300 0,24700
O3 0,33800 0,25600 0,07300
O4 0 0 0
Рис. 1. Структура кристала гідроксиапатиту.
Кластер
Са9(РО4)6
О
Р
Са
ОН
Елементарна
комірка
Са9(РО4)6·Са(ОН)2
53
Присутність кластерів фосфату кальцію розміром від 7 до 10 Å виявлено в
синтетичній рідині організму (СРО) [3] методом динамічного світлорозсіяння (ФКС) та
атомної силової мікроскопії (АСМ). ФКС зарекомендував себе, як досить чутливий
метод для визначення частинок розміром в межах 1 нм. Автори [3] припускають, що
кристал гідроксиапатиту росте з кроком 8 Å, і це пояснюється кластерною моделлю
росту гідроксиапатиту. Методом АСМ підтверджено, що одиницею росту ГАП є кластер
Са9(РО4)6, а не звичайні іони, і ріст кристала відбувається через накопичення цих
кластерів [4]. Встановлено, що кластери Са9(РО4)6 існують в структурі домішок
гідроксиапатиту – октакальцію фосфату та аморфного фосфату кальцію (АФК). Таким
чином, ці кластери можуть бути також структурними одиницями АФК, оскільки вони
забезпечують необхідний зв'язок для переходу фаз АФК → ГАП. Агрегати АФК
змінюються на ГАП за умов впорядкованості їх внутрішньої структури [4]. Нами
зроблено припущення, що ріст кристалу буде обумовлений подальшою взаємодією
утвореного кристала ГАП з Са(ОН)2 і кластерами Са9(РО4)6.
Система диференційних рівнянь кінетики росту кристала гідроксиапатиту
Математична модель динаміки росту кристала в усій повноті свого дослідження
повинна включати диференційні рівняння матеріального, енергетичного балансу та
кінетики [5, 6]. В даній роботі ми обмежились диференційними рівняннями кінетики, що
дають вичерпну відповідь на питання впливу температури та часу синтезу на розмір
наночастинок гідроксиапатиту.
При моделюванні цього процесу, по-перше, було розглянуто ідеальну рівноважну
кінетику, тобто такі процеси, що протікають без порушення статистично рівноважного
розподілу енергії за ступенями свободи. З статистичної термодинаміки відомо, що
реагуючі молекули володіють різною енергією обертальних, поступальних і внутрішніх
ступенів свободи. Ріст кристалів відбувається внаслідок взаємодії молекул з енергією,
рівною чи такою, що перевищує деяке граничне значення – енергії активації. При цьому
доля молекул, здатних до перетворення, повинна в ході реакції зменшуватись. По-друге,
реакція протікає при постійному об’ємі та при відсутності впливу дифузії (створені
умови ідеального перемішування) – реактор періодичної дії.
При дослідженні кінетики в реакторі періодичної дії отримують криві залежності
зміни кількості реагентів перетворення від часу, які називаються кінетичними кривими.
Задача полягає в тому, щоб на основі запропонованої схеми росту кристала ГАП
отримати рівняння, що описує кінетичні криві. Таке рівняння дозволяє обрахувати
кількість речовини в заданий момент часу і при цьому обов’язково повинно містити
початкову концентрацію речовини та постійну, при заданих умовах досліду, величину,
що характерна для даного перетворення (константу швидкості) [7]. Особливість
протікання хімічних перетворень полягає в тому, що утворення продуктів реакції в
помітних кількостях в різних реакціях відбувається за істотно різний час.
Оскільки в ході перетворення об’єм реакційного простору не змінюється, вираз
для швидкості перетворення матиме вигляд
...ba p
b
p
a
a CkC
dt
dC
R == .,
де −ba p
b
p
a CC , концентрації реагентів a, b відповідно; k – константа швидкості хімічної
реакції; pv – стехіометричний коефіцієнт.
Виходячи з цього виразу, сформуємо систему диференційних рівнянь кінетики
росту кристала гідроксиапатиту; дослідження будуть проводитись щодо кристала, який
формуватиметься з 3÷40 кластерів Са9(РО4)6:
54
40
40
)2(
108
8
86
6
65
5
53
3
20
3
53
1
20
3
53
)(
...............................................
...............................................
2,223
,2,3332
kl
kl
iklkli
klі
klkl
kl
klkl
kl
klkl
kl
klkl
kl
j
klmklkl
kl
j
klmklkl
ОНСа
R
dt
dC
RR
dt
dC
RR
dt
dC
RR
dt
dC
RR
dt
dC
RR
dt
dC
jmRRR
dt
dC
jmRRR
dt
dC
=
−=
−=
−=
−=
−=
=−−−=
=−−−=
+
=
=
∑
∑
де і = 3,5,2n, n = 3,...,20, m – кількість кластерів Са9(РО4)6 в кристалі; j = 3,4,5,...,20; Сklі –
концентрація кристалів з і-ю кількістю в ньому кластерів Са9(РО4)6, моль/м
3; ССа(ОН)2 –
концентрація сполуки Са(ОН)2, моль/м
3; Rklm – вираз швидкості росту кристала з m
кількістю в ньому кластерів. Тоді
.
................................................
;
...............................................
;
;
;
;
;
3
Са(ОН)
1
38
2
14040
3
Са(ОН)
1
)2(
2
1
3
Са(ОН)
1
8
2
11010
3
Са(ОН)
1
6
2
188
3
Са(ОН)
1
5
1
166
3
)Са(ОН
1
3
2
155
3
Са(ОН)
3
133
2
2
2
2
2
2
2
CСCkR
CСCkR
CСCkR
CСCkR
CСCkR
CСCkR
CCkR
klklklkl
іklklklіkli
klklklkl
klklklkl
klklklkl
klklklkl
klklkl
=
=
=
=
=
=
=
−
де kkli – константа швидкості росту кристала з і-ю кількістю в ньому кластерів.
В системі диференційних рівнянь у відповідності до методики синтезу [8] задана
початкова кількість речовини (ν) Са(ОН)2 та кластера Са9(РО4)6 становить по 100 моль
при сталому об’ємі реакційної суміші 1 м3. Для цього було прийнято припущення, що всі
реагенти, які були введені, прореагували і за деякий час утворили продукти реакції у
заданій вище кількості. Розв’язок системи диференціальних рівнянь з даними
початковими умовами, що відповідають поставленій задачі, було отримано за
допомогою відомого методу чисельного розв’язку звичайних диференційних рівнянь –
Рунге–Кутта четвертого порядку [9].
Таким чином, отримана система з двадцяти диференційних рівнянь для
дослідження кінетики росту кристалів гідроксиапатиту. Розв’язання такої системи
рівнянь не викликає ускладнень при наявності інформації про величину її параметрів.
(1)
55
Параметри системи диференційних рівнянь
Для визначення термодинамічних параметрів та формування залежностей їх
зміни від температури використаний багатоцільовий пакет квантовохімічної програми
МОРАС [10]. В основу її покладено розв’язок рівняння Шредінгера напівемпіричними
методами квантової хімії в рамках методу молекулярних орбіталей [11].
До параметрів системи диференційних рівнянь кінетики росту кристалу
гідроксиапатиту (1) відносяться константи швидкості росту кристалу з і-ю кількістю
кластерів у ньому. В припущенні, що активований комплекс знаходиться в
термодинамічній рівновазі з реагуючими речовинами, можна отримати такий вигляд
константи швидкості
( ) ( )RTGhTkk ∗−= ∆χ exp1 ,
(2)
де χ – трансмісійний коефіцієнт, приймаємо рівним 1; k1 – стала Больцмана, яка має
значення 1,381·10-38
Дж/K; h – постійна Планка, яка має значення 6,626·10-34 Дж·с;
** STHG ∆∆∆ −=∗ – функція активації Гіббса; ..
*
реагпер НHH −=∆ – ентальпія активації,
Дж/моль; ..
*
реагпер SSS −=∆ – ентропія активації, Дж/(моль·К); .реагН , .перH – ентальпія
утворення реагента та перехідного стану відповідно; .реагS , .перS – ентропія утворення
реагента та перехідного стану відповідно.
Зазвичай, в напівемпіричних методах використовується валентне наближення і
при розв’язку рівняння Шредінгера частина інтегралів опускається чи замінюється на
емпіричні параметри, що відображається на зниженні достовірності та точності
розрахунків. Але у порівнянні з неемпіричними методами, які мають більшу
достовірність результатів розрахунку, напівемпіричні методи мають істотну перевагу в
швидкості їх отримання. Проведення розрахунків неемпіричними методами таких
складних систем, як ГАП, ускладнюється обмеженими можливостями сучасної
обчислювальної техніки по забезпеченню необхідного об’єму пам’яті та швидкодії [12].
Отримання перерахованих параметрів передбачає процедуру мінімізації енергії
реагуючих частинок, які в нашому випадку складаються з кластерів. Мінімальна модель
кристала ГАП – Са3(ОН)6·3Са9(РО4)6, утворюється внаслідок взаємодії кластера
Са3(ОН)6 з трьома кластерами Са9(РО4)6 [2, 4, 13]. Як показує квантовохімічний
розрахунок, проведений методом РМ6, кластер Са3(ОН)6 має структуру з найменшим
значенням повної енергії (рис. 2, а), яка не аналогічна кристалічній структурі
портландіту [14], а рівноважна геометрія кластера Са9(РО4)6 [4, 15] представлена на
рис. 2 б.
а
б
Рис. 2. Оптимізована структура кластерів: а – Са3(ОН)6; б – Са9(РО4)6.
Щоб переконатися в можливості заміни неемпіричних методів розрахунку
напівемпіричними, проведено порівняння геометрії оптимізованих конфігурацій
56
структур гідроксиапатиту, отриманих за допомогою програм МОРАС і РС GAMESS
[16]. На рис. 3 приведені оптимізовані конфігурації структур кластера Са9(РО4)6
напівемпіричним методом РМ6 і неемпіричним методом Хартрі-Фока-Рутаана з
базисним набором 6-31G*. Як видно, структури цього кластеру, отримані
напівемпіричним і неемпіричним методами, подібні між собою та із структурою,
представленою на рис. 1.
а б
Рис. 3. Оптимізована структура кластера Са9(РО4)6: а – напівемпіричний метод РМ6; б –
неемпіричний метод 6-31G*.
Провівши аналогічну процедуру оптимізації геометрії для продукту реакції
взаємодії трьох кластерів напівемпіричним і неемпіричним методами (рис. 4) та
порівнявши їх молекулярні структури, переконуємось в допустимості використання
напівемпіричного методу РМ6 при дослідженні таких систем.
а б
Рис. 4. Оптимізована структура мінімальної моделі кристала гідроксиапатиту: а –
напівемпіричний метод РМ6; б – неемпіричний метод 6-31G*.
На рис. 5 приведений один з імовірних механізмів початку росту кристала ГАП.
Але це можливо тільки у випадку наявності деякого надлишку енергії (енергії активації)
у порівнянні з мінімальною енергією реагуючих частинок (кристал, кластер). Для
визначення цієї енергії необхідно знайти стабільний стан реагента, який буде
відповідати мінімальному значенню енергії, та конфігурації перехідного стану цього
процесу на поверхні потенціальної енергії. Складність знаходження їх полягає в тому,
що поверхня потенціальної енергії багатоатомної моделі має, як правило, окрім
глобального мінімуму, також і велику кількість локальних мінімумів. Ітераційні
57
алгоритми оптимізації, які використовуються в сучасних квантовохімічних програмах,
не дозволяють напевне відшукати глобальний мінімум систем з великою кількістю
атомів. Знайдений в результаті оптимізації екстремум може бути глобальним, а може і
не бути ним.
Тому в даній роботі для пошуку глобального мінімуму проводилась оптимізація з
різних початкових положень як кластера Са9(РО4)6, так і кластера Са3(ОН)6, стосовно до
вже утвореної частинки (кристала) з наступним порівнянням значень енергії в точках
різних локальних екстремумів. Умовою закінчення пошуку мінімуму енергії в
градієнтних методах оптимізації є мінімальна величина градієнта функції відгуку.
Вихідні реагенти Продукт реакції
Рис. 5. Механізм початку росту кристала гідроксиапатиту.
Після знаходження конфігурації молекулярної моделі, що відповідає реагенту і
продукту, знаходиться конфігурація перехідного стану реагуючої системи. Існують
спеціальні алгоритми для пошуку перехідного стану, наприклад алгоритм слідування
власному вектору чи алгоритм синхронного переходу [1]. Такі алгоритми передбачають
знаходження сідлової точки. Але чим складніша система, тим важче знайти однозначний
розв’язок, що пов’язано з наявністю великої кількості стійких станів на потенціальній
поверхні. Знайдена в результаті критична точка може бути, як сідловою точкою так і
точкою локального мінімуму. Слід також зазначити, що тривалість часу пошуку сідлової
точки значно перевищує тривалість часу пошуку глобального мінімуму. У зв’язку з цим
пошук перехідного стану здійснювався методом релаксаційного сканування
потенціальної поверхні. Такий алгоритм, звісно, вже заздалегідь вносить деяку похибку,
але така похибка буде незначною через багатовимірність потенціальної поверхні, де
різниця між локальним мінімумом і максимумом невелика.
Для перевірки вірності такого підходу був проведений розрахунок знаходження
геометрії перехідного стану кристала з 5 кластерів напівемпіричним методом NLLSQ
програми МОРАС [10]. За стартову точку молекулярної структури взято конфігурацію
системи, отриманої за описаним вище алгоритмом. Результат, засвідчив незначну
відмінність їх ентальпій активації, близько 5 кДж/моль.
Розрахунки матриці Гессе за допомогою програми МОРАС дали змогу отримати
термодинамічні параметри: ентальпії та ентропії реагентів та перехідних станів реакції
при різних температурах.
Са3(ОН)6
58
Слід зауважити, що знаходження рівноважних структур кристалів з вибраною
нами кількістю кластерів (40 кластерів), навіть за допомогою неемпіричного методу, є
неможливе. Тому розрахунки кристалів із нарощуванням кількості кластерів
проводились до тих пір, поки різниця між поточним та кількома попередніми
значеннями термодинамічного параметра кристала з 40 кластерами, які отримуються за
еквівалентною лінією (побудованою відповідно за різної кількості значень параметрів
згідно з методом середньоквадратичного відхилення) залежності ∆Н*=f(і) не ставала
меншою за деяку допустиму величину. Наприклад (рис.6), виходячи з наявності 7 точок,
отриманню яких передувало 7 розрахунків кристала з 5, 6, 8, 10, 12, 14, 16 кластерами у
ньому, побудована пунктиром еквівалента лінія залежності ∆Н* = f(і) при температурі
353 К. При врахуванні кожної з наступних останніх трьох точок, для побудови
еквівалентної лінії, різниця значень по ній ентальпії активації для кристала з 40
кластерами не перевищувала 5 кДж/моль.
Рис. 6. Розрахункові залежності зміни ∆Н* в функції кількості кластерів в кристалі.
Для підстановки в вираз константи швидкості (2) ентальпії активації у
відповідності з графічними даними, представленими на рис.6, виведено поліноми залеж-
ності ∆Н*=f(Т) для кристалів з усіма кількостями кластерів (3–40). Для прикладу, на
рис. 7 по розрахунковим даним (рис. 6) побудовані характеристики залежності ∆Н* = f(Т)
(пунктирна лінія) та виведені поліноми лінійної апроксимації, згідно з якими побудовані
залежності у вигляді суцільної лінії для кристалів з 10 та 40 кластерами відповідно. За
таким же алгоритмом знайдено поліноми залежності ∆S* = f(Т) для кристалів з такою ж
кількістю кластерів.
а б
Рис.7. Розрахункові залежності зміни ∆Н* в функції температури Т кристала з: а –10
кластерів; б – 40 кластерів.
59
Нижче наведені вирази розрахунку константи швидкості для кристалів з 10 та 40
кластерами, у кожному з виразів присутні поліноми залежності ∆Н*=f(Т) та ∆S*=f(Т), які
були отримані по вище наведеному алгоритму.
kkl10 = (1,381-23·T/6,626-34) ·exp(-((- 0,21·T + 1,72)×1·103-T·(- 0,53·T + 70))/(8,314·T));
kkl40 = (1,381-23·T/6,626-34) ·exp(-((- 1,12·T + 4,22)×1·103-T·(- 0,53·T + 40))/(8,314·T)).
Аналогічно одержані вирази константи швидкості для кристалів з іншою кількістю
кластерів.
Результати та обговорення
Таким чином, сформувавши систему диференційних рівнянь кінетики росту
кристалу та отримавши на основі розрахунків методами квантової хімії параметри до
неї, а саме константи швидкості росту кристала було проведено чисельний експеримент
синтезу нанодисперсного гідроксиапатиту.
Розроблена математична модель кінетики росту кристала гідроксиапатиту на
даному етапі дає змогу дослідити вплив температури синтезу і концентрацію реагентів в
розчині на розміри частинок. На рис. 8 – 11 приведені результати розрахунку впливу
температури синтезу на динаміку зміни кількості кристалів гідроксиапатиту за
розмірами в 1 м3 реакційної суміші. Розміри кристала визначаються кількістю кластерів
у ньому. Так, позначення кривих згідно з (1) відповідають такому алгоритму
(3Кр = 3Р+О, 5Кр = 5Р+2О, (2*n)Кр = (2*n)Р+nO, n = 3,…20),
де 3Кр, 5Кр,(2*n)Кр, n = 3,…20 – нумерація кристалів з відповідною кількістю у ньому
кластерів Са9(РО4)6 та Са3(ОН)6; Р – кластер Са9(РО4)6; О – кластер Са3(ОН)6.
а
б
Рис. 8. Кінетичні криві зміни кількості: а – вихідних реагентів; б – кристалів
гідроксиапатиту за розмірами при температурі синтезу 278 К.
а
б
Рис. 9. Кінетичні криві зміни кількості: а – вихідних реагентів; б – кристалів
гідроксиапатиту за розмірами при температурі синтезу 298 К.
60
а
б
Рис. 10. Кінетичні криві зміни кількості: а – вихідних реагентів; б – кристалів
гідроксиапатиту за розмірами при температурі синтезу 318 К.
а
б
Рис. 11. Кінетичні криві зміни кількості: а – вихідних реагентів; б – кристалів
гідроксиапатиту за розмірами при температурі синтезу 353 К.
При температурі синтезу 278 К, в розчині переважатимуть кристали
гідроксиапатиту з 6, 5 і 10 кластерами (з розмірами в ~5 – 8 нм) з кількістю речовини
приблизно 1,75; 1,0; 0,7 моль відповідно. Концентрація кристалів з більшою кількістю
кластерів є дуже незначною. Вже при температурі синтезу 353 К значно переважатиме
кількість (1,4 моль) кристалів з 40 кластерами (з розмірами ~32 нм). Отримані
результати свідчать про те, що підвищення температури синтезу веде до зростання
середньої енергії вихідної системи та зростання числа частинок, з необхідним
надлишком енергії, і, внаслідок цього, до підвищення швидкості реакції. І головне до
збільшення розмірів кристалів. Такі результати досить непогано співвідносяться з
експериментальними даними [17], а також з результатами робіт [18, 19], які свідчать, що
із зростанням температури синтезу зростають розміри кристалів в розчині (рис.12).
Достовірність отриманих результатів і висновків забезпечується:
Рис. 12. Залежність розмірів кристалів від
температури синтезу [18].
61
- відповідністю використаних фізичних наближень і математичних моделей
квантової хімії поставленій задачі;
- коректністю квантовохімічних обчислювальних методів;
- задовільним співпадінням розрахункових значень, отриманих на основі напів-
емпіричних розрахунків геометричної структури і розрахунків чисельними методами
диференційних рівнянь, з наявними експериментальними даними.
Висновки
Таким чином, розроблено математичну модель кінетики росту кристалів
гідроксиапатиту, яка дає змогу проаналізувати вплив зміни кількості вихідних реагентів
і температури синтезу на розміри частинок. Отримані результати добре корелюють з
результатами фізичного експерименту. Виходячи з розрахунку, оптимальна температура
та час синтезу утворення нанодисперсного гідроксиапатиту з розмірами 3÷5 нм
становить 278 K та 1÷4×104
с, відповідно.
Розроблена математична модель та отримані на її основі результати можуть в
подальшому бути використані при розробці розширених моделей на основі чисельного
методу розв’язку задач прикладної фізики та хімії – методу кінцевих елементів. Такі
моделі одночасно враховуватимуть процеси, які пов’язані з масоперенесенням, транс-
портуванням енергії і з хімічною кінетикою, та можуть ефективно використовуватись,
при створенні та дослідженні кераміки і композиційних матеріалів з біополімерною ма-
трицею, а також засобів доставки біологічно активних речовин і терапії на клітинному
рівні.
Література
1. Соловьев М.Е., Соловьев М.М. Компьютерная химия. – M.: СОЛОН-Пресс. – 2005. –
536с.
2. Fiorentini V., Methfessel M. Extracting convergent surface energies from slab
calculations // J. Phys. Condens. Mater. –1996. – V.8, N 36. – Р.6525–6529.
3. Onuma K., Ito A. Cluster growth model for hydroxyapatite // Chem. Mater. – 1998. –
V.10. – Р.3346–3351.
4. Onuma K. Recent reseach on pseudobiologycal hydroxyapatite crystal growth and phase
transition mechanisms // Progr. in Cryst. Growth Character. Mater. – 2006. – V.52. –
Р.223–245
5. Кубасов А.А. Химическая кинетика и катализ. – Москва: Изд-во Моск. ун-та. –
2004. – 302с.
6. Эмануэль Н.М., Кнорре Д.Г. Курс химической кинетики. – Москва: Высш. шк. –
1984. – 463 с.
7. Эткинс П. Физическая химия. – Москва: Мир. – 1980. – 1145с.
8. Hayek E. Newsely pentacalcium monohydroxyorthoposphate (hydroxyapatite) // Inorg.
Synth. – 1963. – V.7. – P. 63–65
9. Ануфриев И.Е. Самоучитель Matlab 5.3/6.x. – СПб.: БХВ-Петербург. – 2003. – 736 с.
10. Аминова Р.М. Расчеты электронного строения и свойств молекул
полуэмпирическими методами квантовой химии. – Казань – 1997. – 71с.
11. Кобычев В.Б. Квантовая химия на ПК: Компьютерное моделирование молекулярных
систем. – Иркутск: Иркут. гос. ун-т. – 2006. – 87 с.
12. Кобзев Г.И. Применение неэмпирических и полуэмпирических методов в кваново-
химических расчетах. – Оренбург: Оренбург. гос. ун-т. – 2004. – 150с.
13. W. Zhu, P. Wu. Surface energetics of hydroxyapatite: a DFT study // Chem. Phys. Lett. –
2004 – V. 396 – Iss.1–3. – P. 38–42
14. Desgranges L., Grebille D., Calvarin G., Chevrier G., N. Floquet, Niepce J.-C. Hydrogen
thermal motion in calcium hydroxide: Ca(OH)2 // Acta Cryst. – 1993. – V.49. – P. 812-817.
62
15. W. J. McCarthy, D. M. A. Smith, L. Adamowicz, H. Saint-Martin, I. Ortega-Blake / An аb
initio Study of the Isomerization of Mg− and Ca−Pyrophosphates // J. Am. Chem. Soc., –
1998. – V.120. – N.24. – Р. 6113–6120.
16. Schmidt M.W., Baldridge K.K., Boatz J.A., Elbert S.T., Gordon M.S., Jensen J.H., Koseki
S., Matsunaga N., Nguyen K.A., Su S.J., Windus T.L., Dupuis M., Montgomery J.A.
General atomic and molecular electronic-structure system: Review // J. Comput. Chem. –
1993. – V. 14. – N 11. – P.1347–1363.
17. Головань А.П., Туров В.В., Барвинченко В.М., Мищенко В.М., Горбик П.П.,
Шевченко Ю.Б. Наноструктурированные композиты на основе белков костной ткани,
высокодисперсного кремнезема и гидроксиапатита // Хімія, фізика та технологія
поверхні. – 2007. – Вип. 13. – С. 309.
18. Бакунова Н.В., Баринова С.М., Шворева Л.И. Влияние температуры синтеза на
размер наночастиц гидроксиапатита // Рос. нанотехнологии. – 2007. – Т. 2. – №9–
10. – С.102–105.
19. Pang Y.X., Bao X. Influence of temperature, ripening time and calcination on the
morphology and crystallinity of hydroxyapatite nanoparticles // J.Eurор. Ceram. Soc. –
2003. – V.23. – P.1697–1704.
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ КИНЕТИКИ РОСТА
КРИСТАЛЛА ГИДРОКСИАПАТИТА
А.П. Головань1, И.В. Головань2, Е.Н. Демяненко
1
1
Институт химии поверхности им. А.А. Чуйко Национальной академии наук Украины,
ул. Генерала Наумова, 17, Киев, 03164, Украина
2
Институт электродинамики, Национальной академии наук Украины,
пр. Победы 56, 03680 Киев, Украина
Разработана математическая модель кинетики роста кристаллов гидроксиапатита,
которая дает возможность проанализировать влияние изменения количества вещества
исходных реагентов и температуры синтеза на размеры частиц. С помощью
квантовохимических расчетов получены массивы зависимостей термодинамических
параметров от температуры и выведены их полиномы. Проведен численный эксперимент и
получены характеристики динамики изменения количества частиц по размерам для различных
температур синтеза.
MATHEMATICAL MODELLING OF THE HYDROXYAPATITE
CRYSTAL GROWTH KINETICS
А.P. Golovan1, I.V. Golovan2, E.M. Demianenko1
1Chuiko Institute of Surface Chemistry, National Academy of Sciences of Ukraine,
17 General Naumov Str. Kyiv, 03164, Ukraine
2Institute of an electrodynamics, National Academy of Sciences of Ukraine,
Peremogy Prosp. 56, 03680 Kyiv, Ukraine
A mathematical model of hydroxyapatite crystal growth kinetics which enables to analyse
influence the change of initial reagents substance quantity and synthesis temperatures on particle sizes
is developed. Arrays of dependences of thermodynamical parameters on the temperature are obtained
with quantum-chemical calculations and their polynoms are deduced. Numerical experiment is carried
out and characteristics of dynamics changes of quantity particles in the sizes for various synthesis
temperatures are obtained.
|