О математических особенностях связанных моделей теплои массопереноса в твердых средах

О математических особенностях связанных моделей тепло- и массопереноса в твердых средах

С.Н. Сорокова, А.Г. Князева

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

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

On mathematical features of coupled models of heat and mass transfer in solids

S.N. Sorokova and A.G. Knyazeva Institute of Strength Physics and Materials Science SB RAS, 634021, Tomsk, Russia

We study coupled systems of equations of heat and mass transfer for solids by the example of models of metallothermic transformations and diffusion in a binary system. The coupled systems are shown to be systems of variable type equations; the type of systems is determined by the domain of model parameters. A comparative analysis of the models leads to the conclusion that heat and mass transfer processes in solids are not similar in the general case.

1. Введение

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

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

теплопереноса и деформирования в твердой фазе при распространении металлотермических превращений

[1]. Вторая модель учитывает взаимовлияние диффузии и деформирования при диффузии в бинарной системе

[2]. Обе модели анализируются в приближении малых деформаций и содержат «подобные» малые параметры — время релаксации потока тепла и время релаксации потока массы; а также малые параметры, характеризующие различие скоростей процессов переноса и распространения механических возмущений. В обеих моделях выполняются условия для компонент тензора деформаций £ у и напряжений Сту:

£11 Ф 0, £22 = £33,

£кк = £11 + £22 + £33 ~ £11 = £, £12 = £31 = £23 = 0, ст11 Ф 0, ст22 = СТ33 ф о, ст12 = ст23 = ст31 = 0,

т.е. реализуется случай «одноосной» деформации.

2. Распространение фронта превращения в металлотермической системе

Математическая постановка задачи о распространении фронта превращения в металлотермической системе следует из [1, 3] и включает:

© Сорокова С.Н., Князева А.Г, 2005

- уравнение теплопроводности с химическим источником, связанное с деформациями:

ЭТ Э Эе Эу

сер— = ——Jт — ЗКатТ— + Q-е Эt Эх т т Эt ^

Эt’

- уравнение химическои кинетики для суммарной реакции А ^ В

Э- = К Ф( у) ехР |^—^, Ф( У) = 1 — у,

- уравнение движения или баланса импульса Э 2 и да,, рэ^=1^

Модель замыкают соотношения между напряжениями и деформациями, следующие из модели Максвелла:

11 + —ап = (Л + 2ц) — К—— + — К (е — м),

Эt

Эа22 , Ц к

дt

Эе ^ Эм ц — К , дt Эt к

Эе — К Эм + ц

Эt Эt к

22 + —а22 — Л^ — К+ £ К (е — м),

где w = 3 [аТ (Т - Т0) + (ав -а А)у], и уравнение для потока тепла в виде обобщенного закона Фурье:

Jт — Лт

ЭТ ЭJ7

— и

дх дг

Граничные условия зависят от конкретной ситуации и, разумеется, не оказывают влияния на тип систем уравнений.

В (1)-(5) приняты следующие обозначения: г—время; х — пространственная координата; гг — время релаксации теплового потока; р — плотность; с£ — теплоемкость; АТ — коэффициент теплопроводности; JT — вектор плотности теплового потока; Qг — суммарное тепловыделение (поглощение тепла) в результате реакции; аТ — линейный коэффициент теплового расширения; у — доля суммарного продукта реакции или степень превращения; кг — предэкспоненциальный множитель; Еа — энергия активации суммарной реакции; R — универсальная газовая постоянная; аА, ав — коэффициенты «концентрационного» расширения реагента и продукта; А, ц — коэффициенты Ламе; К — изотермический модуль всестороннего сжатия (К = А + +2 ц/3); к — коэффициент динамической вязкости.

Введение переменной V = дм/дг позволит представить систему (1)-(5) в виде системы шести дифференциальных уравнений перового порядка по времени. В этом случае вместо (1) и (3) имеем:

Эт _ ЭJт Эи Эу

Се р— —-----------------ЗКатт------+ Ог —,

е Э* Эх т Эх Э*

Эу Эа11

р^- —

Эt Эх

д£ = дv дг дх

Для уменьшения числа параметров системы уравнений (2), (4)-(8) запишем ее в безразмерных переменных:

(т — Т»)Еа ^ „ х е

9= ' а, т —-, £ ——, е — —,

ЯТ»2 и х» е»

а

а»

где

3»

х» —д/к^» , кт — ,

СеР

V

V»

t» —— ехр

ко

ЯТ„

СеРКт»

Qr Еа к0

ехр

ЯТ„

3» —

2ц + Л

Лт ят»2

, V» — -

е» х»

х» Еа t»

представляют собой удобные для задачи масштабы, выбранные из физических соображений [3]. В результате получим систему уравнений вида:

_Э0 _ Эт

ЭЯт

Эт

— ш(0 + р)^- +

ЭК 1 Эу

э^+7эГ ’

1 Э9 ят

тг Э£ тг ’

Эу л ч

— У(1 — у) ехр

Эт

Р0 +1

1 Э5

^ Э^ ’

зг

Эт

Э5 5 —

Эт + ту

— Эг — Э0_ ^Эу 2_ Э^ Эт у Эт ту

Эе — ЭГ

эГ—"Э^ ’

h е — (9+90)— — у У

где у —

серЯТ»2

яп

2 — р кТ .

Ц 7 2 К

ту — ^ t»; h —

к

(Т» — То)Еа ЯТ»2

t» 2ц + Л t-»

а в — а а .

т , "1 ’ ^ I ; 90

2ц + Л атQо/cеР

параметры модели (первые четыре —

малые параметры в соответствии с их физическим смыслом). Если ю = 0, имеем несвязанную задачу, если а = 0, имеем квазистатическую задачу.

Для определения типа системы дифференциальных уравнений, линеаризуем уравнение теплопроводности в окрестности температуры продуктов 0в, которая находится в процессе решения задачи [3].

Полученная в итоге система уравнений представима в матричном виде:

^+в ^ = г,

дт д£,

0

где U = (0, qT, у, V, s, е) — вектор искомых величин; В — матрица коэффициентов, свойства которой определяют тип системы уравнений; F — вектор правых частей.

Уравнение |B - А1| = 0 позволяет найти собственные значения матрицы В.

Если а = 0, тг = 0 система уравнений сводится к параболической задаче, т.е. к решению уравнения теплопроводности совместно с уравнением кинетики.

Если а = 0, тг Ф 0, то решение задачи сводится к решению уравнения теплопроводности гиперболического

типа: А, 2 = ± -------1------- и А3 4 5 6 = 0.

12 у Тг [1 + ю(0в +Р)] 3’4’5’6

Если а Ф 0, тг = 0, получаем: А12 = ± 1/ а и А3 456 = = 0. В этом случае система уравнений имеет парабологиперболический тип.

В общем случае имеем А, 2 = 0, остальные собственные значения следуют из решения биквадратного уравнения:

тга2А4 - [а2 + тг

т.е. тип системы уравнений зависит не от всех параметров модели. Характерная область изменения параметров, влияющих на тип системы, такова: тг = 0 10, а = 0 1.5, ю = 02, в = 0.3. Это уравнение не имеет комплексных корней. Если а Ф 0, тг Ф 0, то во всей области изменения параметров система уравнений (2), (4)-(8) имеет гиперболический тип.

Если система уравнений несвязанная, то уравнения разных типов решаются независимо.

3. Изотермическая диффузия в бинарной системе

Одномерная математическая модель диффузии в бинарной системе в полубесконечном слое (например, при насыщении его с поверхности легирующим элементом или при ультразвуковом воздействии) следует из [2, 4] и включает уравнение диффузии

ЭС1 — —

Эt Эх ’

уравнение движения (3) и уравнение для потока массы в виде обобщенного закона Фика, учитывающего конечность времени релаксации потока массы и влияние напряжений на массоперенос:

J — —д

ЭС1+В

Эх р

С

Эа

кк

Эх

-т

Э£

г Эt

Компоненты тензоров напряжений и деформаций связаны соотношениями

Сту = 2ц£г>- + 8у (А£кк - 3К(а1 - а2)С1),

С2 + С1 = 1.

В (10)-(12) С1 — концентрация легирующего эле-

а ^

мента; В0 = ——1—1 — коэффициент переноса под RT *

действием напряжений; Д = D-[g11 — коэффициент

диффузии элемента в бинарной системе; £п — термо-

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

вий имеем £„ = £ =------ и

11 дх

кк = СТ11 + 22 + СТ33 =

11 22

33

— (3Л + 2ц)е — 9 К (а1 — а 2)С1.

Эи

Тогда после введения переменной V = — система

дг

уравнений (3), (10), (11) примет вид:

дv . д£ дС1

^—= (А + 2ц)-— 3К (а1 -а 2)—^-, дг дх дх

ЭС1 — — Э£ Эt Эх ’

ЭJ ж

tr — г Эt

: — і —

д*ё , 9К(а1 — а2) ВС

д1 ёи +----------------ВоС1

р

ЭС1

Эх

3КВоС1 Эе

р Эх

Эе — Эу Эt Эх

В безразмерных переменных

с — с гї — J г, — е V — у £ — х т — t С — С1, Q — т , Є — , Г — тг , £— , т— ,

3» е» V» х» и

где х» =

д/д1"ї»;

а» е» х»

; е» —---; а» — 3Ка.; V» —--------------------.

Л + 2ц 1 t»

х»

3» — —, имеем

и

а2 ЭГ Эе а ЭС ЭС ЭQ Эе ЭГ а "Э7—э£— "Эр "Эт — —'ЭТ эГ —"Э^’

Э0 /-> ЭС 2 Эе

т^3~ — —й — (Ян + юС

Эт Э£ Э£

где тг — 1»; а2 —

р х»2 а1 — а 2

-; А — —--------------------------------; ю —

а

9Ка1 от1 ; h2 —

К

— параметры модели. Безраз-

рЯТ ' Л + 2ц

мерное время релаксации тг, малый параметр а и коэффициент связанности имеют аналоги в предыдущей модели.

Для определения типа системы уравнений массоуп-ругости линеаризуем уравнение для потока массы в окрестности некоторого характерного значения концентрации С, при достижении которого вклад в суммарный массоперенос «бародиффузии» становится сравним со вкладом концентрационной диффузии и которое находится в ходе решения задачи [5, 6]. Далее поступаем аналогично. Представив полученную систему уравнений в векторной форме:

ЭU . ЭU _

---+ А— — F,

Эт Э£

,

Рис. 1. A = -0.5, Tr = 0.01, область: ш = 0.001

Рис. 2. A = -0.5, а = 0.1, область: ш = 0.001

Рис. 3. A = -0.5, ш = 0.1, область: тг = 0.005

В обоих частных случаях в указанной области изменения параметров система уравнений имеет парабологиперболический тип.

В практически интересных частных случаях [7] значения параметра А могут далеко выходить за пределы исследованной области. Тогда при некотором соотношении параметров модели и величины концентрации С возможно появление комплексных корней А34 и проявление эллиптичности связанной системы уравнений массоупругости. Это возможно, если ~п + юС -- Аюh2 < 0.

В общем случае (а Ф 0, тг Ф 0) можно выделить малую область параметров, в которой система уравнений имеет эллиптический тип. На рис. 1-3 ей соответствует область между кривыми 1-1,2-2, 3-3; при увеличении тг область эллиптичности увеличивается и сдвигается вверх по а. При увеличении коэффициента связанности область эллиптичности в параметрах (а, С) или (тг, С) системы увеличивается.

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

Работа выполнена при финансовой поддержке РФФИ, грант № 03-01-00074.

где U = {V, С, е, 0} — вектор искомых функций; А — матрица коэффициентов; Р — вектор правых частей, вычислим определитель | A -А1| = 0. В результате имеем биквадратное уравнение для определения собственных значений:

тга2А4 - [тг + (^~11 + юС)а2]А2 +

+

Уравнение исследовано в широкой области изменения параметров: А = -0.5, 0, 0.5, h = 0.9, ю = 0.001^0.1, а = 0.00001^0.1, тг = 0.00001^0.1. Было принято

<11 =1

В частном случае тг = 0, а Ф 0 получаем: А12 = 0

и X 3 4 = ±.

~ii + шС — АшН2 \ (~ii + шС)а2

В частном случае а = 0, тr Ф 0 получаем: X12 = 0

g~11 + шС — АшН2

Литература

1. Князева А.Г., Дюкарев Е.А. Модель распространения стационарного фронта превращения в вязкоупругой среде // ФГВ. - 2000. -Т. 36. - № 4. - С. 41-51.

2. Князева А.Г. Перекрестные эффекты в твердых средах с диффузией

// ПМТФ. - 2003. - Т. 44. - № 3. - С. 85-99.

3. Князева А.Г., Сорокова С.Н. Асимптотический анализ задачи о рас-

пространении безгазового горения в вязкоупругой среде // Физ. мезомех. - 2004. - Т. 7. - Спец. выпуск. - Ч. 1. - С. 62-65.

4. Князева А.Г. Диффузия и реология в локально-равновесной термо-

динамике // Сб. научных трудов «Математическое моделирование систем и процессов». - Пермь: ПГТУ, 2005. - № 13. - С. 45-60.

5. Князева А.Г., Тимохин А.М. Стационарные режимы распростране-

ния простейшей твердофазной реакции, лимитируемой диффузией // Сб. «Химическая физика процессов горения и взрыва»: XII Симпозиум по горению и взрыву. - Черноголовка: ИХФ РАН, 2000.- Ч. III. - С. 93-95.

6. Князева А.Г. Режимы развития из начального зародыша твердофазной реакции, лимитируемой диффузией // ФГВ. - 1996. -Т. 32. - № 4. - С. 72-76.

7. Князева А.Г., ДонскаяЯ.Г. Диффузионно-деформационная модель

развития зародыша продукта твердофазной реакции // ФГВ. -1997. - Т. 33. - № 2. - С. 52-68.

т