Численное исследование диффузионных процессов в имплантатах с многослойными биоактивными покрытиями при их взаимодействии с модельной биологической жидкостью Текст научной статьи по специальности «Физика»

Численное исследование диффузионных процессов в имплантатах с многослойными биоактивными покрытиями при их взаимодействии с модельной биологической жидкостью

Н.Н. Назаренко, А.Г. Князева, И.А. Хлусов1, А.В. Карлов1

Институт физики прочности и материаловедения СО РАН, Томск, 634021, Россия 1 Центр ортопедии и медицинского материаловедения СО РАМН, Томск, 634050, Россия

Предложена модель диффузионного взаимодействия имплантата с многослойным покрытием и физиологического раствора. Эта модель представляет собой сопряженную задачу многокомпонентной диффузии. Принято условие малости коэффициентов диффузии в твердой фазе по сравнению с жидкостью, что позволило пренебречь распределением концентраций элементов в физиологическом растворе. Проанализированы частные варианты модели. Определены поля концентраций в твердой фазе и среднее содержание концентраций в растворе в различные моменты времени при варьировании параметров модели. В реальной ситуации варьирование коэффициентов диффузии достигается за счет изменения состава покрытий.

Numerical study of diffusion processes in implants with many-layered bioactive coatings at interaction with a model biological liquid

N.N. Nazarenko, A.G. Knyazeva, I.A. Khlusov, and A.V. Karlov

In the work a mathematical model of diffusion interactions between an implant with a multilayered coating and physiological solution is proposed. The model is the conjugate problem of multicomponent diffusion. The condition of the smallness of the solid phase diffusion coefficient in comparison with the diffusion coefficient in a liquid phase allows us to neglect the distribution of concentrations in the physiological solution. Special variants were analyzed in the model. Concentration fields in the solid phase and their middle concentrations in the solution were found for different time moments, varying the problem parameters. In real situation, a change of the diffusion coefficient is reached by a change of the phase and chemical component of coatings.

1. Введение

Кальций-фосфатные материалы и покрытия имплантатов, относящиеся к классу биоактивных материалов, находят все более широкое распространение в травматологии и ортопедии благодаря высокой способности к интеграции с костной тканью. Они делятся, в свою очередь, на поверхностно активные (гидроксиапатит, бифазная керамика (гидроксиапатит + трикальций-фос-фат), стеклокерамика) и резорбируемые (рассасывающиеся) материалы на основе сульфатов и фосфатов кальция. Концепция биоактивности рассматривает высвобождение в окружающую среду (в результате биодеградации искусственного материала) биологически активных ионов (ионов кальция и фосфора, ионов легирующих примесей и т.д.), которые влияют на остеоинтеграцию имплантатов.

Биодеградацию имплантата можно моделировать in vivo и in vitro за счет физико-химического и биохимического растворения материала в агрессивных жидкостях и/или посредством биорезорбции, осуществляемой клеточными системами организма (макрофаги, остеоклас-

ты). Однако чувствительность анализов in vitro выше, чем in vivo. При этом концепция биоактивности применима как для объемных керамических имплантатов, так и покрытий на металлах. Биодеградация объемных кальций-фосфатных материалов хорошо известна. Растворимость многослойных кальций-фосфатных покрытий на металле более сложна, изучена в меньшей степени. В то же время, построение уравнений движения ионов на границе покрытие/жидкая среда позволяет прогнозировать поведение имплантата в костной ткани (расшатывание или врастание), и следовательно, успех лечебного процесса.

2. Формулировка задачи

Математическая модель диффузионного взаимодействия имплантата на основе титана с многослойным покрытием рассматривается в следующей постановке (рис. 1). Полагаем, что процесс можно считать одномерным, поскольку толщины диффузионных слоев и покрытий много меньше размеров (диаметра) исследуемого образца; пренебрегаем образованием химических

© Назаренко Н.Н., Князева А.Г., Хлусов И.А., Карлов А.В., 2004

Рис. 1. Иллюстрация к общей постановке задачи

соединений, выделением новых фаз и влиянием внутренних напряжений на процесс диффузии.

Начальное распределение основных элементов (Л, Са, Ка, С1) считаем известным. Для массовых концентраций в твердой фазе вводим обозначения С1, С2, С3, С4 соответственно. Все остальные элементы в твердой фазе обозначим через С5. В жидкой фазе концентрации элементов обозначим через Аг. Баланс масс записывается следующим образом:

Е с = 1, Е а, = 1. (1)

г-1 г-1

Каждое покрытие представляет собой многофазную и многокомпонентную систему, может иметь пористую структуру, что требует привлечения дополнительных соображений. В данной задаче пренебрежем перекрестными диффузионными потоками, но учтем, что коэффициенты диффузии каждого элемента зависят от концентраций всех элементов. Тогда полная постановка задачи примет вид:

^ =-D1L'q (А.) ^ + Я,, х < 0, ді дх 1 К дх 1 (2)

х = -кь: - DlLIq ^ = 0, ох (3)

х = 0, DLq ^ = Da ^, А = С,, ох ох (4)

дС, д дС, 1Т -ъ1* С) их •х 5 0- (5)

х = к,: Du = 0, ох (6)

і = 0: Сі = С0,Аі = А0. (7)

границах раздела имеем:

(Рі )ь ^ = (Di )р ^, С1Ъ = Ср. ох ох (8)

В (1)-(8) индекс i относится к номеру элемента, j — к номеру слоя; индекс L относится к величинам слева от границы раздела, а индекс Р — к величинам справа. Зависимости коэффициентов диффузии от концентрации соответствуют теории Вагнера для многокомпонентных сплавов. Источник (или сток) аг каждого элемента может быть связан с разбавлением или осаждением элемента на стенках. В последнем случае = К гСг, где Кг — коэффициент осаждения, i = 1, 2. Учитывая, что коэффициенты диффузии в жидкой фазе на 3-5 порядков выше, чем в твердой фазе, пространственным распределением концентраций в жидкос-

ти пренебрежем. Тогда, проинтегрировав уравнения диффузии для жидкой фазы (2) в пределах от -hL до 0 с учетом условий на границе раздела фаз (4) и условия (3), придем к задаче твердофазной диффузии (5)-(8) с условием:

х = 0: - D^

С

дх

= к

с

ді

(9)

В задаче требуется исследовать поведение концентраций элементов в объеме жидкой фазы во времени, а также их перераспределение в покрытии для различных параметров, входящих в модель.

Принимая некоторые вполне разумные упрощения, можно сделать аналитические оценки, полезные для обработки данных эксперимента и оценки некоторых констант, входящих в модель. Частные варианты модели позволяют исследовать подробно влияние некоторых факторов на процесс диффузионного перераспределения элементов, не прибегая к анализу полной модели.

3. Примеры частных моделей

3.1. О возможности оценки коэффициентов диффузии

Диффузионная задача вида

- = Dk

ді к

х = 0: - Dk

дх 2

дСк

(10)

да

х = — D

дх $Ск с дх

= к^, ак (0, і) = Ск (0, і), (11)

ді

= 0,

(12)

I = 0: Ск(х, 0) = 0, ак(0) = а0 (13)

дает распределение (до разбавления) концентраций тех элементов, которые в начальный момент времени содержались только в жидкости с концентрацией а°. Приближение «бесконечной» твердой фазы для малого времени диффузии (неделя для диффузии в твердой фазе — малое время) и большая толщина твердого слоя не будут играть особой роли, что можно показать специальным образом. Задача легко решается операционным методом и дает распределение натрия или хлора в жидкости:

ак = ак ехр

Dt

ег&

(14)

где ег&(z) = —-г= Г е у dy. л1п І

При условии і >> кь/Dk, что для реальных значений коэффициентов диффузии выполняется всегда, функция ег&(і) имеет асимптотическое представление, используя которое вместо последней формулы, найдем

ак = а° Д . (15)

Л/лDkі

Следовательно, обрабатывая данные эксперимента по концентрации этих элементов в жидкой фазе (в раст-

Рис. 2. Иллюстрация к возможному определению коэффициента диффузии

воре) в переменных С1к''^‘ (по оси ординат) и (по

аХ 41

оси абсцисс), можем оценить коэффициенты диффузии этих элементов в твердой фазе по углу наклона полученной «прямой» (рис. 2).

Обрабатывая данные по содержанию элементов в жидкости для каждого из твердых материалов отдельно, определим эффективные коэффициенты диффузии нужных элементов для каждого слоя многослойного покрытия. Более того, если проводить специальный эксперимент с твердыми материалами, содержащими разный процентный состав каких-либо элементов, на основе простейшей модели (10)-(13) принципиально возможно определить зависимость эффективных коэффициентов диффузии от концентраций этих элементов.

Аналогичная модель, имеющая аналитическое решение, имеет место и для тех элементов, которые в начальный момент времени находятся только в покрытии.

3.2. Диффузия в пористом материале

Рассмотрим диффузию элементов (натрия и хлора) из физиологического раствора в пористый материал (рис. 3), который считаем однослойным. Математическая постановка задачи диффузии в такой системе может быть представлена в форме:

дС]_ = D д С1 дt

(С1 - С2),

ЭС,

дt

= D^

х = 0: D1

* _ а

дх2~ ~

+ ^ (С1 - С2),

Эх2 Л

ЭС1 , ЭС1

■ = к

Эх

дt ’

(16)

(17)

(18)

ЭС1

Эх

= к

дС2

дt ’

дС1 дС2

х ^^: —1 = 0,—- = 0,

Эх Эх

(19)

(20)

t = 0: С1 = С2 = С0, х = 0,

С1 = С2 = 0, х > 0, (21)

где С1 — концентрация элемента в матрице; С2 — его концентрация в порах; D1 — коэффициент диффузии в матрице; D2 — коэффициент диффузии в порах; DL — коэффициент диффузии в жидкости (физиологическом растворе), причем DL >> D1, DL >> D2, D = D1/D2 << << 1, а — коэффициент обмена между фазами; Л — пористость; С0 — начальное содержание элемента в жидкости. Эта задача также имеет аналитическое решение, но оно не удобно в использовании.

Численное исследование задачи показывает, что увеличение коэффициента обмена между фазами, очевидно, приводит к увеличению содержания диффундирующего элемента в матрице. Пример показан на рис. 4, где введены обозначения

£ £

т = - ’ £ = -

г, а = а £

л/5=-

Рис. 3. Иллюстрация к постановке задачи диффузии в пористый материал

£*52’ ’

где и» — характерное для задачи время. Чем меньше 5, тем быстрее убывает концентрация элемента в жидкой фазе

4. Численное исследование полной модели

Для численного решения задачи (5)-(9) использованы переменные £, = х/к0, где к0 —толщина верхнего покрытия, равная 200 мкм, и т = £/£», и» — период разбавления, коэффициенты диффузии обезразмерены по коэффициенту диффузии титана в слое 4. В задаче имеются два основных параметра 51 = £5»/(характерное расстояние, которое «пробегает» титан в результате самодиффузии за время £*) и Ь = к1^/к0 , условием применимости модели будет 51 << 1, что в реальном диапазоне изменения параметров всегда выполняется.

Проведено подробное численное исследование задачи в широкой области изменения параметров. Анализ результатов показал, что зависимостями многих коэффициентов диффузии от концентраций можно пренебречь, что существенно сокращает число параметров, которые необходимо найти из эксперимента. Этот результат связан с малыми концентрациями Ка и С1 всюду, Т в областях 1, 2 и Са в областях 3, 4. Качественно результаты получаются одинаковыми, как для коэффициентов диффузии, не зависящих от концентрации, так и зависящих от них.

При учете разбавления раствора в два раза в моменты времени х = 1, 2, 3,..., то есть каждую неделю, для этих моментов времени в точке х = 0 задавались условия:

Рис. 4. Распределение концентрации в матрице и в порах образца. ц = 0.3, С0 = 0.1, В = 0.001, = 100. Сплошная линия соответствует

расчету с а = 0, пунктир — а = - 0.01; символы — а = — 0.001

Рис. 5. Распределение концентрации титана в системе в различные моменты времени (^ < т2 <^з < ^4) и зависимость концентрации титана в жидком растворе от времени для различных значений коэффициента DTІ в верхнем покрытии < вТ2 < вТ3)

Рис. 6. Распределение концентрации кальция в системе в различные моменты времени (^ <т2 <^3 <Т4) и зависимость концентрации кальция в жидком растворе от времени для различных значений коэффициента 5Са в верхнем покрытии (5^^* < 5(?а < 5(С3а))

Сі = Сі/2, С2 = С^2,

С3 = 0.014 + С3/2, С4 = 0.0215 + С4/2, где СТІ, ССа, СКа, СС1 — решение, полученное к данному моменту времени. Примеры расчетов представлены далее.

На рис. 5 показано распределение концентрации титана в твердой фазе к моменту времени % = 6 (т.е. через шесть недель) и изменение его содержания в жидкости для различных коэффициентов диффузии титана в верхнем слое. Видим, что если все коэффициенты диффузии одинаковы и равны единице (кривая 1), несмотря на разбавление концентрация титана в физиологическом растворе растет. При увеличении же коэффициента диффузии титана в 5 раз (кривая 2) устанавливается некоторый квазистационарный режим. При увеличении коэффициента диффузии титана в 10 раз (кривая 3) концентрация титана в растворе сначала возрастает, а затем уменьшается.

Увеличение коэффициента диффузии кальция в верхнем слое (рис. 6) в пять раз приводит вначале к возрастанию его содержания в растворе (кривая 2), а

при увеличении в 10 раз к уменьшению (кривая 3), что связано с исчерпанием кальция в покрытии. На концентрациях других элементов это не сказывается. Кривые 1 соответствуют расчету с одинаковыми коэффициентами диффузии.

Варьирование параметров модели показало, что осаждение элементов на стенках не является основным фактором, способствующим уменьшению их концентрации в растворе.

Варьирование коэффициентов диффузии в реальной ситуации достигается изменением фазового и химического состава покрытий.

5. Заключение

Моделирование процессов диффузии элементов покрытия в биологических жидкостях позволяет предсказывать поведение биоактивных имплантатов в организме, следовательно, уже на технологическом этапе формировать изделия с требуемыми биомедицинскими свойствами. Подобные разработки имеют теоретический, практический и экономический эффект.