Методы исследования диффузионных процессов в твердых телах на основе волновой модели диффузии Текст научной статьи по специальности «Математика»

В.А. Сиренек

УДК 666.11:541.64

МЕТОДЫ ИССЛЕДОВАНИЯ ДИФФУЗИОННЫХ ПРОЦЕССОВ В ТВЕРДЫХ ТЕЛАХ НА ОСНОВЕ ВОЛНОВОЙ МОДЕЛИ ДИФФУЗИИ

(Санкт-Петербургский государственный технологический институт) e-mail: wasirenek@gmail.com

Разработаны математические методы (аналитические и численные, детерминированные и вероятностные) исследования свойств волновой модели диффузии в применении и расчету релаксационных процессов массопереноса.

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

Существующие (как правило, эмпирические) способы учета релаксационного характера массопереноса (кажущегося, с позиций уравнения Фика, эффекта "запаздывания" диффузионных потоков) нельзя считать удовлетворительными, т.к. они не приводят к теоретическому обоснованию и, тем более, к точному расчету основных стадий диффузии. Автором разработана единая методология математического моделирования релаксационных процессов массопереноса в гетерогенных системах "твердое тело - реагент", основанная на использовании гиперболических уравнений волновых (релаксационных) моделей [1]. Наиболее простое из них имеет форму телеграфного уравнения. Гиперболические уравнения, в отличие от параболических с постоянными коэффициентами (в том числе от уравнения диффузии Фика), учитывают конечность скорости распространения возмущений концентрации. Их решения позволяют адекватно описывать дифференциальные характеристики массопереноса, т.е. профили концентрации целевого компонента с резко очерченным фронтом, проявляющим себя на начальной (релаксационной) стадии процесса. (Решения уравнения Фика адекватно описывают стадию развитой диффузии с "размытым" фронтом профиля концентрации). Из экспериментов по изучению кинетики диффузии в системах "твердое тело - реагент" чаще всего известны интегральные характеристики массопереноса (количество переносимого вещества, ширина диффузионной зоны), способы практического использования которых для расчета дифференциальных характеристик не достаточно разработаны.

Гиперболическое уравнение (волновая модель) диффузии:

тс; + 4= ^ (1)

является композицией выражения для потока и уравнения неразрывности:

д = -0*с'х -т4, (2)

с;+д'х = 0, (3)

где q - поток, О* - эффективный коэффициент массопереноса, т - время концентрационной релаксации. Использованы асимптотические формы уравнения (1) - волновое (при «<<т) и параболическое (при «>>т) уравнения (далее - индексы "в", "п"), определяющие кинетические законы лимитирующих стадий диффузии.

Многие приложения теории волновой модели диффузии приводят к задаче о продвижении концентрационного фронта в полуограниченном теле. Задача расчета начальной стадии диффузионного извлечения целевого компонента из образца на основе волновой модели и ее решение рассмотрены в безразмерной форме:

2"тт + 1'т= ; 7(X,0) = 1, 1'т(X,0) = 0, Х>0;

г(0, Т)=0, 7>0, (4)

7(х,т )4 - X р' '^У^/2) л - е~х/2 ПрИ X < Т; 1 прих > Л (5) I 2 ^ Л2 - X2 ]

где Ъ = (с - Сгр)/(Сн - Сгр); Сн, Ср, с(х, «) - концентрации извлекаемого компонента: начальная (в толще образца), на границе с реагентом, текущая; Сн>Сгр; Х=х/(О*т)1/2; Т = «/т.

Нулевая скорость изменения концентрации компонента в начальный момент времени обусловлена твердой фазой. Профиль 2(Х, Т) отражает существование в момент Т невозмущенной области (Х>Т) и диффузионной зоны (Х<Т) (рис. 1). Скачок концентрации на фронте профиля уменьшается во времени по закону ехр(-Т/2). Практически при Т>10 профиль 2(Х, Т) с "размытым" фронтом может быть адекватно аппроксимирован решением уравнения Фика.

Для изучения кинетики развития диффузионной зоны в твердых телах использована средняя эффективная ширина диффузионной зоны Н(Т); при этом Л(«)=Н(«/т)(О*т)1/2 - размерная

величина. Для задачи (4) получены выражения:

H(T) = (1 - Z(X, T))dX = JT e-S!2I0 (S / 2)dS =

= Te T/2 (I0 (T/2)+1, (T/ 2)),

(6)

где I0, I - модифицированные функции Бесселя.

Рис. 2. Средняя эффективная ширина диффузионной зоны

H(VT) для различных моделей диффузии Fig. 2. Average effective wight of diffusion zone H(VT) for different models of diffusion

Рис. 1. Профиль концентрации Z(X,T) (график решения (5) задачи (4))

Fig. 1. Profile of Z(X,T) concentration (curve of solution of (5) for task (4))

Первое из (6) удобно при анализе мгновенной Жмгн, а второе - средней Жср скорости роста H. Значение T*=t*/T=4/rc«1.3 определяет время смены асимптот функции H(T) - законов лимитирующих стадий диффузии: HB(T)=T и НП(Т)=(2/%1/2)Т1/2 (точка M на рис. 2); при этом H*=H(T)™1, h*=h(t*)«(D*x) 1/2. Функция H(T1/2) при Т>Т* выходит на участок слабой нелинейности (квазистационарный режим), где может быть достаточно хорошо аппроксимирована отрезком (CB ) прямой (AB), выходящей из точки (A), сдвинутой вправо от начала координат. Это позволяет объяснить наблюдаемый на практике в координатах (h, t1/2) эффект "запаздывания" диффузионных потоков. Временной области эффективного применения волновой модели диффузии (1), т.е. значениям t < 10т, соответствует интервал значений: h(t)<h(10T)»3.5(D*T)1/2.

Закон Максвелла - Каттанео (2), положенный в основу уравнения (1), будем трактовать как следствие закона ёд/& = - (д - д0)/т для скорости возвращения (релаксации) некоторой характеристики д системы, выведенной из состояния равновесия, к исходному значению д0 [2]. В качестве д рассмотрим поток вещества, а в качестве д0 - его стационарное значение в виде 1-го закона Фика. При этом параметр т определяет время релаксации потока д , т.е. время, за которое отклонение д от д0 уменьшается в е раз. В развитие этого подхода изучены релаксационные свойства различных характеристик скорости роста Н(Т). Вычислены периоды релаксации (Тр=р/т) величин Ж2мгн(Т) и Жср(Т), отражающие свойства волновой модели диффузии. Так, для Ж 2мгн значение Тр«1.2 служит оценкой Т*, а для Жср значение Тр«9 является оценкой временной области применимости модели.

Исследована связь гиперболических уравнений массопереноса и их решений со случайными имитационными процессами, т.е. процессами, лежащими в основе вывода этих уравнений [3]. Это предоставляет возможность не только решать основанные на волновых моделях массопереноса диффузионные задачи методом Монте-Карло, но и изучать вероятностную природу моделируемых физических явлений. Проведено вероятностное исследование кинетических характеристик волновой модели диффузии. Получен вероятностный аналог выражения (5):

2(Х,Т)=?{\^а (Т )|<Х}={Т )(Х) ,

UT) = \T (-1)Na (s)ds, « =1/2 ,

(7)

где ^а(Т) - смещение "блуждающей" со скоростью ^ =1 "частицы" за время Т; Е\^(Т)\(Х) - функция

распределения величины |^а(Т)\, а - интенсивность (частота) перемены "частицей" направления движения. Ширину диффузионной зоны субъективно "привязывают" к различным уровням приведенной концентрации ^=0.4 0.8. С учетом (7) такой величиной для 2и=р будет р - квантиль Ир функции Е, т.е. ^а(Т)\(Ир) = Р (рис. 1). Исследованы два вида средних значений Ир, не зависящих от выбора р. Использована теорема о том, что случайная величина и ее Р - квантиль совпадают по распределению вероятностей, если р - случайная величина, равномерно распределенная на (0,1). Величина Е{Ир} определяет среднюю эффективную ширину диффузионной зоны И, а (Е{И2р})1/2 -среднеквадратическое смещение "частицы" ("диффузионный путь" Идиф):

Нд

h (t ) = £ и в ^ = б{яв }= ешг)|}= 0 (1 - % (Г, иф(?) = ^ И^Р^ Б{Ир } =7 Б{^(Г)} = 72(T -1 + exp(-T)).

(8)

(9)

Расчеты по (7) и (8) выполнены методом Монте-Карло при 105 реализациях Ъ,а(Т). Расхождение результатов расчета по (5) и (7), а также по (6) и (8) для Т< 5 и Х< 5 - менее 1%, для Т< 5 и Х<0.5 - менее 0.1%. Заметим, что Идиф^*)»И*& «(Б*т)1/2, Лдиф(10т)«4.2(Б*т)1/2. При О>т для Б*=Б из формулы (9) следует соотношение Лдиф(0«(Б*01/2, соответствующее уравнению Фика; при К<т получаем Лдиф(0«(Б*/т)1/2^ Характер зависимостей Н(() и Лдиф(?) для лимитирующих стадий диффузии проявляется и в эмпирических оценках процесса формирования диффузионной зоны с учетом граничной химической кинетики [4].

Вероятностная интерпретация характеристик Z(X, Т) и Н(Т), согласие их теоретических и опытных значений (в наших расчетах диффузии в металлах и стеклах в [1] средняя относительная ошибка 5 <10%) дают основание для использования уравнения (1) и его стохастического аналога при описании процесса массопереноса на микроскопическом уровне. Элементарный цикл модели "случайного блуждания", соответствующей уравнению (1), состоит из движения "частицы" по прямой со скоростью ^ от одного момента изменения направления движения до следующего. Продолжительность цикла и соответствующее ему перемещение Sw - одинаково показательно распределенные случайные величины, средние значения которых приняты за микропараметры: <и>=1/а, <Sw>=w/a. Макропараметрами служат коэффициенты т и Б* уравнения (1): т=<^/2,

Б*=<«У/2^.

В идеальных моделях перемещения кинетических единиц (атомов, групп атомов, ионов или молекул) в неупорядоченных структурах твердого вещества для ансамбля частиц за микропараметры принимаются - среднее время пребывания их в метастабильном состоянии (в химической связи со структурой вещества) <т'> и среднее расстояние между локализациями таких состояний <А>. Временем перескока из одного метаста-бильного состояния в другое пренебрегают. За макропараметр принимается коэффициент диффузии из уравнения Фика Б=<А>2/<т'>.

Полного соответствия между обеими схемами перемещения частиц на микроуровне нет. Введены соотношения пропорциональности между

их микропараметрами: <^>=^^4^, <Х>=к2<«м>. Критерием соответствия схем на макроуровне будем считать сопоставимость макропараметров Б и Б* при условии, когда гиперболическое уравнение можно аппроксимировать параболическим, при этом ^=2к22. По поводу соотношения т«<т'> следует заметить, что отождествление времени релаксации процесса миграции частиц со средним временем их "оседлой жизни" проводилось давно [5]. Коэффициент Б уравнения Фика связывает оба микропараметра и не позволяет идентифицировать их по отдельности. Волновая модель диффузии позволяет устранить этот недостаток.

В [6] разработан двухэтапный метод расчета коэффициентов уравнения (1) по данным о выходе целевого компонента из образца. На первом этапе вычисляются оценки. Предложена аппроксимация Н(Т) в классе оценок и (г ;ц) типа среднего гармонического НВ(Т) и

Нп(Т): (1/И)ц = (1/Ив)ц + (1/Ип)ц - со средней относительной ошибкой 5 (Г; ц). Для значения Ц=2: 5(Г;2) <, (И(ТГ;2) - кривая 4 на рис. 2). Для ц=2.55: 5(Г;2.55) <1%; при этом для Т<Т*: 5(Г ;2.55) < 0.1%.

С учетом аналитических выражений для НВ(Т) и НП(Т) получено:

0С) =

(сн - С]

гр

(к )12[гц 2 + (т/ к )ц/ц

Б* =

к

Л2п2/ Ц

1. \2/ц

Ь ^ Ц, к = я/4.

(сн Сгр) а т= к2 ^ , к = я/4, (10)

где и (Сн - Срр - количество переносимого

вещества. Величины а и Ь определяются как коэффициенты уравнения регрессии, следующего из (10): /=ац + Ь, / = ^=гц/2. Формулы для а и

Ь получены на основе N опытных данных аналитически методом наименьших квадратов (МНК) для абсолютной (А) и для относительной (в) задаваемых погрешностей уд измерения Q; нелинейность координат относительно 0 учтена в МНК с "весами":

1

N

тт 52(а, Ь) ; 52 =-1--V1

аЬ (N - 2) ^

I=1

А - аЧ1- Ь . У ( / )

у( /) =

5/

•Уе

- погрешность /.

На втором этапе значения Б* и т уточняются. Оценка т служит начальным приближением к корню (11), полученному при решении МНК задачи о наилучшей аппроксимации экспериментальных значений Qэ(ti) модельной Q(t):

2

e;> W ) - о«, >!• q«, ) -^f;,(;; "Д' -».

;-1 /2т) + /2т)]

где

Q(t;) = V D 4; (11)

4 = (t^VT)exp(-;/2т)[/0 (/,/2т) + ^ (/,/2т)];

4 оэ(;)/ZNi

А 2 =1А .

В области 1>1* (при выходе данных процесса в координатах (д^) на квазистационарный режим) правомерна линейная аппроксимация: Q(1) = а1112 + Ь1. Уравнению Фика соответствует Q(1) = а2112. Расчетные формулы для а1, Ь1, а2 получены для тех же видов погрешностей измерения Q, что для а и Ь из (10). Коэффициент диффузии В для этих двух моделей рассчитывается по формуле:

В = ка12/(еи - сгр)2, а1=^а, (12 )

а - угол наклона прямой при аппроксимации данных в координатах (д). Формула (12) получена из условия: г(х) = 1и , где г(х) = сгГ(х/2)- автомодельное решение уравнения Фика; ~ = х . Каждому значению 2и соответствует свое к. Для 2и=0.5 к =1; для 2и=0.8: к =1/я. Средней эффективной ширине диффузионной зоны для уравнения Фика соответствует Д^0.6 и к =л/4 [7].

Ранее в [7] в качестве меры толщины гипотетического слоя, ответственного за эффект

"запаздывания" диффузионных потоков, исполь-

0 1 /2 2 2 зовалась величина к ~(В1з) , где 1з=Ь1 /а1 - мнимое время "запаздывания" ((0А)2, рис. 2). Автором в качестве такой характеристики принято к*«(В*т)1/2; при этом к*>к0.

Для различных видов задаваемых погрешностей измерения Q (или к) (абсолютной и относительной) и различных моделей диффузии в зависимости от стадии процесса (волновая модель, уравнение квазистационарного режима, уравнение Фика) разработан алгоритм и программа для расчета систем "твердое тело - реагент" [8]. Примеры расчета реакционной диффузии в металлах и процессов выщелачивания стекла с учетом релаксационного эффекта массопереноса приведены в [1].

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

ЛИТЕРАТУРА

1. Таганов И.Н., Сиренек В.А. Волновая диффузия. СПб.: НИИХ. 2000. 209 с.;

Taganov I.N., Sirenek V.A. Wave diffusion. SpB: NIIKh. 2000. 209 p. (in Russian).

2. Сиренек В.А., Щеголев В.В. Новый справочник химика и технолога. СПб.: НПО Профессионал. 2004. С. 296-302; Sirenek V.A., Shchegolev V.V. New handbook of chemist and technologist. SpB: NPO Professional. 2004. P. 296-302 (in Russian).

3. Сиренек В.А. // Физика и химия стекла, 2002. Т. 28. №3. С. 221-229;

Sirenek V.A. // Phyzika i khimiya stekla. 2002. V. 28. N 3. P. 221-229 (in Russian).

4. Гегузин Я.Е. Диффузионная зона. М.: Наука 1979. 343 с.; Geguzin Ya.E. Diffusion zone. M.: Nauka. 1979. 343 p. (in Russian).

5. Френкель Я.И. Введение в теорию металлов. М.: Физ-матгиз. 1958. 368 с;

Frenkel Ya.I. Introduction to metal theory. M.: Fizmatgiz. 1958. 368 p. (in Russian).

6. Сиренек В.А. // Физика и химия стекла. 2003. Т. 29. № 4. С. 507-519;

Sirenek V.A. // Phyzika i khimiya stekla. 2003. V. 29. N 4. P. 507-519 (in Russian).

7. Белюстин А.А., Шульц М.М // Физика и хим. стекла. 1983. Т. 9. № 3. С. 3-27;

Belyustin A.A., Shults M.M. // Phyzika i khimiya stekla. 1983. V. 9 N 3. P. 3-27 (in Russian).

8. Сиренек В.А., Соков В.М. Программа определения параметров диффуз. процессов на границе твердое тело-реагент. ОФАП. Свид. о регистр. №8801 от 05.04. 07; Sirenek V.A., Sokov V.M. Programm for determination of diffudion process parameters on solid-reagent interface. OFAP. Registration certificate N 8801 from 05.04. 07 (in Russian).

Кафедра системного анализа