CC BY
70
Алгоритм численного решения задач неизотермической диффузии, встречающихся в процессах поверхностной обработки
Н.В. Букрина, А.Г. Князева
Институт физики прочности и материаловедения СО РАН, Томск, 634021, Россия
Предложен алгоритм численного решения задачи неизотермической диффузии, учитывающий различие пространственных и временных масштабов процессов переноса в твердой фазе. Применимость алгоритма продемонстрирована на примере задачи остывания нанесенного покрытия, содержащего легирующие элементы. Показано, что учет зависимости теплопроводности от температуры сказывается на перераспределении элементов в диффузионной зоне. Модель имеет приложение к процессам поверхностной обработки материалов с использованием высокоэнергетических источников.
Numerical solution algorithm for non-isothermal diffusion problems in surface treatment processes
N.V. Bukrina and A.G. Knyazeva Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634021, Russia
The paper proposes a numerical solution algorithm for a nonisothermal diffusion problem, which accounts for the difference in spatial and temporal scales of transfer in the solid phase. The algorithm validity is illustrated by the example of a problem on cooling of a deposited coating with alloying elements. It is shown that an account of the temperature dependence of thermal conductivity affects element redistribution in the diffusion zone. The model can be applied in processes of surface treatment of materials using high-energy
1. Введение
Покрытия на металлические изделия наносят для повышения механических, антикоррозионных и других свойств поверхности этих изделий или для сохранения заданных свойств в условиях воздействия внешней среды.
При термических методах нанесения покрытий происходит локальный разогрев присадочного материала и подложки, иногда до расплавления, способствующий интенсивному молекулярному и конвективному переносу компонентов. Аналогичные процессы наблюдаются при высокотемпературной поверхностной обработке, где локальное расплавление может быть вызвано экзотермическими химическими превращениями. В любом случае в формировании свойств материалов важную роль играют диффузионные процессы, приводя к перераспределению легирующих элементов.
При численном решении задач, моделирующих подобные технологические процессы, возникает ряд проб-
лем, обусловленных различием характерных скоростей и пространственных масштабов различных физикохимических явлений в конденсированной фазе. Это приводит к необходимости как построения моделей, учитывающих особенности технологических условий, так и к разработке специальных численных алгоритмов [1]. Например, двумерная математическая модель формирования диффузионной зоны в неизотермических условиях при условии диффузии одного легирующего элемента представлена в [2]. Насыщение поверхности металлов несколькими элементами позволяет обеспечить такие технологические свойства, которые нельзя получить при поверхностном легировании одним элементом. При этом возникает еще одна проблема — определение парциальных коэффициентов диффузии и их зависимости от температуры и концентраций.
Цель настоящей работы заключается в разработке специального численного алгоритма для решения зада-
© Букрина Н.В., Князева А.Г., 2006
Рис. 1. Иллюстрация к постановке задачи
чи неизотермической диффузии с учетом перекрестных диффузионных потоков.
2. Математическая постановка задачи
Математическая постановка задачи об остывании системы «покрытие - подложка» включает в себя уравнение теплопроводности и уравнения диффузии легирующих элементов из покрытия толщиной к1 в основу (рис. 1):
ср-
дТ
дt да дt да
дt
1 __
2 __
где
31 - -Вц
32 --В 21
■—Л—
дх дх ’ д31 дх ’ д3 2 дх ’
да1
(1)
(2)
дх
да1
В
12
- В
да2 дх ’ да2
дх 22 дх
есть диффузионные потоки легирующих элементов; В у , I = 1, 2, — парциальные диффузионные коэффициенты; а { — концентрации легирующих элементов; Т — температура; с, р, Л — теплоемкость, плотность и коэффициент теплопроводности, зависящие от температуры и концентрации диффундирующих веществ.
Примем, что толщина слоя Н1 много меньше, чем толщина подложки, что позволяет не рассматривать распределение температуры и концентрации в этом слое. Это условие корректно, если выполняются неравенства:
Нх г * В1, г, у -1,2,
к < у г * к ,
где к-Л^(с1р1); г * — характерное время процесса; В*- — коэффициенты диффузии элементов в покрытии.
Тогда в качестве граничного условия для температуры при х = 0 будет справедливо граничное условие в виде дифференциального уравнения:
х-0: Л— -А1с1р1 — + ае(Т4 -Те4),
дх
дг
(3)
где слагаемое ае(Т4 - Те4) описывает теплопотери по закону Стефана-Больцмана; Те — температура окружающей среды.
Для концентраций на границе имеем аналогичные условия:
х - 0: к
да1
дг
-- 315 к1
да2
дг
- - 3 2
(4)
На бесконечном удалении от нагретой области источники и стоки тепла и массы отсутствуют:
дТ
х : Л—-0, дх
^ - 0, ^ - 0.
дх дх
(5)
(6)
В начальный момент времени распределение концентраций и температуры в покрытии однородно:
г -0:
Т - Тх,
а^ — аю, а 2 — а
20
х - 0,
Т - Т
0,
х > 0.
а1 - 0, а2 - 0
Концентрация третьего элемента следует из соотношения а з -1 - а1 - а 2.
Учитываем, что теплофизические свойства и коэффициенты диффузии зависят от температуры. Зависимости теплофизических свойств от температуры аппроксимируются полиномами. Коэффициенты диффузии в общем случае зависят от концентраций по степенному закону, а от температуры — по закону Аррениуса:
Ву - (ВУ + + В?1а 2 + В|а1а2 +
+ ВУа1 + В5а 22)ехр
У
■Е-1
RT ''
/
Л-Л 0 +Л1Т + Л 2 Т2 + Л 3Т3 + Л4 Т4
(7)
(8)
Теплоемкости веществ резко возрастают в окрестности температур плавления, причем покрытие и подложка в общем случае имеют различные температуры плавления ТрЬд, ТрЬ. Например, для покрытия имеем
[(ср^,1 + ^рs,1 8(Т - ТрЬ,1), Т < ТрЬ,1,
(ср)1 -
(ср)Щ,
Т > Т
(9)
где 8 — дельта-функция Дирака; индекс s относится к твердой фазе; L — к жидкости. Соотношение, аналогичное (9), имеет место и для подложки.
В расчетах 8-функция Дирака в (9) заменяется 8-образной функцией Ф 0 [3]:
Фл
Ф0 -'
1
ехр
Чл/П
ехр
(Т - ТрЬ)2
(Т - Тръд)
ДЛЯ ОСНОВЫ И покрытия соответственно, где s, S1 — параметры сглаживания, которые подбираются таким образом, чтобы описать характерное плато на температурной кривой в окрестности температуры плавления.
Заметим, что матрица диффузионных коэффициентов должна удовлетворять определенным требованиям. Так, в рассматриваемой трехкомпонентной системе должно выполняться неравенство
D11D 22 — D12 D 21 > 0 следующее из термодинамических ограничений.
Для практических целей интерес представляют трех-компонентны1е системы1 типа Fe+Cr+C, для которые можно найти все параметры [4], требующиеся для расчета.
В задаче требуется найти распределение концентраций и температуры в различные моменты времени в процессе остывания системы.
3. Описание численных алгоритмов
3.1. Задача теплопроводности
В расчетах приняты следующие значения параметров: h1 = 0.05 см; Cs = 0.68 Дж/(г • K), CL = 0.65 Дж/(г • K), рs = 7.61 г/см3, рL = 7.865 г/см3, Lph = 250 Дж/г, Lph1 = = 400 Дж/г, Tph1 = 1 710 K, Tph = 1 810 K.
В простейших предельных случаях задача имеет точное аналитическое решение. Так, при условии постоянства теплофизических свойств и при отсутствии плавления имеем для задачи теплопроводности [5]:
1 л/к
K7T
= T0 + (Tm - T0)exp
X
VK
х exp
Гґ і—л І 1 VK t erf /" К 1 + ^1 I
1 ^ J Кc hi 2J к t \ У
(10)
Но совпадение численного решения с аналитическими в этих предельных случаях — слабое доказательство корректности численного решения нелинейной задачи. Поэтому описанию численного алгоритма уделим особое внимание.
Численное решение задачи теплопроводности не встречает особых затруднений. Решение проводится по неявной консервативной разностной схеме с использованием метода прогонки. Выбор пространственного Ах и временного Аг шагов обусловлен, в основном, необходимостью корректного описания температурной кривой в окрестности температуры плавления.
Так, на рис. 2 представлена зависимость температуры в точке х = 0 от времени.
В расчете без учета плавления кривая Т(г) является гладкой (кривая 1).
Если ЬрЪ Ф 0 , при некотором соотношении параметра сглаживания и шага Аг на кривой Т(г) появляется плато, соответствующее температуре плавления (ср.
Рис. 2. Зависимость температуры поверхности от времени (х = 0); Дх = 0.001; At = 0.0001 с. Кривая 1 — расчет без плавления; кривые 2,3 — расчеты с учетом плавления; результаты для различных значений параметра сглаживания s = Si = 5, 10, 25; результаты совпадают. Кривые 1, 2 — теплопроводность не зависит от температуры; кривая 3 — X = X(T); 4 — Дх = 0.001; At = 0.01 с; s = s1 = 10; X = const
кривые 2 и 4). Плато появляется (при данном наборе физических параметров) как с уменьшением пространственного шага при постоянном s, так и с изменением параметра сглаживания s при постоянном шаге при условии, что на область резкого изменения теплоемкости в окрестности температуры плавления приходится не менее 10 точек. Зависимости объемной теплоемкости С = ср от температуры в окрестности температуры плавления, соответствующие разным значениям параметра сглаживания, показаны на рис. 3.
Требование к шагу по времени можно сформулировать следующим образом: при изменении температуры от Т1 до Т2 за время г2 - г1 в окрестности температуры плавления (рис. 2) должно выполняться неравенство
11 - *2-> 10.
Аг
Эти же рассуждения относятся к величине пространственного шага Ах в области больших градиентов температуры.
Кривые 1, 2 рассчитаны для постоянной теплопроводности X = 0.45 Вт/(см • К • с). Кривая 3 соответствует расчету с теплопроводностью, зависящей от температуры. Зависимость принята в виде
X = 0.72 -1.37-10 -7 Т -
- 6.00-10-7 Т2 + 3.39 • 10-10 Т3.
Плато в окрестности температуры плавления сохраняется, но время остывания при учете зависимости теплопроводности от температуры изменяется существенно. При учете зависимости коэффициента теплопроводности от температуры шаги разностной сетки, приводящие к приемлемому результату, могут оказаться иными.
Рис. 3. Зависимость объемной теплоемкости от температуры: а — Ах = 0.001; At = 0.0001; 5 = — = 5 (1), 10 (2), 25 (3); б — Ах = 0.001; А* = 0.0005; 5 = - = 5 (1), 10 (2), 25 (3), 50 (4)
Представленные результаты подтверждают известное мнение, что в случае нелинейных задач использование абсолютно устойчивых разностных схем не гарантирует получение качественно верного результата при произвольном соотношении шагов разностной сетки.
3.2. Диффузионная задача
При условии независимости коэффициентов диффузии от концентраций диффузионная задача также допускает аналитические оценки [6], из которых следует, что распределение концентраций элементов может быть немонотонным при определенном соотношении парциальных коэффициентов диффузии.
Численно задачу диффузии можно решать различными способами.
Казалось бы, что мелкие шаги по пространству и по времени, обусловленные особенностями решения задачи теплопроводности с плавлением, позволяют использовать для диффузионной задачи простой алгоритм, основанный на явной разностной схеме. На первый взгляд, применение явной разностной схемы для задач диффузии, также предполагающей использование мелкого шага по времени, не внесет дополнительных ограничений. Но для того чтобы на ширину диффузионной зоны приходилось достаточное число точек, нужно, чтобы АхБ г * Б у , что затрудняет решение задач на
единой разностной сетке, поскольку условие устойчивости явной разностной схемы
°шіп^ < 0.5
дх 2
(11)
требует еще более мелких шагов по времени. Это существенно замедляет счет и перегружает память. Поэтому такой вариант численного алгоритма оказывается неприемлемым.
Был разработан специальный алгоритм, использующий особенности процессов теплопроводности и диффузии, а именно учитывающий различие их характерных пространственных масштабов (рис. 4). Решение задач теплопроводности и диффузии осуществлялось
на разных разностных сетках. В диффузионной задаче использована явная разностная схема.
Этот вариант алгоритма также оказывается неудачным, так как не избавляет от ограничения (11). Проблема усугубляется тем, что в неизотермических условиях величина шіп(0^) заранее неизвестна, что затрудняет предварительную оценку.
Поэтому задачу неизотермической диффузии решали с использованием неявных разностных схем для обеих подзадач.
Неявная разностная схема для задачи диффузии, которая решается на своей пространственной сетке т = Аґ/Аґд раз за один шаг по времени Аґ задачи теплопроводности, также может быть составлена различным образом. Использование полностью неявной схемы предполагает применение для полученной системы линейных разностных уравнений метода матричной прогонки. Такой выбор может оказаться неудачным, если коэффициенты диффузии отличаются более чем на порядок. Поэтому было использовано «расщепление» диффузионных уравнений на слагаемые, часть из которых аппроксимировалась разностями с нижнего слоя, т.е. по явной схеме. Это аналогично схемам расщепления по физическим факторам, но отличается от обычных схем расщепления [7]. Такой способ аппроксимации диф-
Рис. 4. Пояснение к согласованию разностных сеток для численного решения задачи неизотермической диффузии
т, к - 0 С1 -0.02 ■ 1 б с2 - 0.16 ■
2000 - *1 «\ 0.01 ■ \\Ч 0 08 '
1000 ■ \ .
\ V ^ -1 — 0.00 ■ . ■ >^,>4 0.00 -
0.00 0.05 0.10
X, см 0.0
И
1.2 -10^ X, см 0.0
1.2-10“4 X, см
Рис. 5. Распределение температуры в образце в последовательные моменты времени 0.001, 0.0025, 0.015 с (а) и концентраций в диффузионной зоне в момент времени 0.001 с (б, в) в процессе остывания материала с покрытием. Штрих-пунктирная кривая — расчет без учета перекрестных диффузионных потоков и зависимости теплопроводности от температуры, X = 0.45 Вт/(см • с • К); штриховая кривая — расчет с учетом перекрестных потоков, X = 0.45 Вт/(см • с • К); сплошная кривая — расчет с учетом перекрестных диффузионных потоков и зависимости
теплопроводности от температуры. Б>1°1 = 1; = 10 1;
= 10-1;
Б2 = 5 • 10-1; Е1
- Е2 = 140 000 Дж/моль. Для задачи теплопроводности Ах =
= 0.001; Аг = 0.00025; для задачи диффузии Ах = 5 • 10 7; Аг = 1.25 • 10
ференциальных уравнений разностными обусловлен еще и тем, что никакие линейные преобразования не позволяют в нелинейной задаче разделить уравнения диффузии подобно тому, как это делается в задачах с постоянными коэффициентами диффузии [8] (см. Приложение).
Пример пространственного распределения концентраций элементов в момент времени 0.001 с и температуры в различные моменты времени показан на рис. 5. Видно, что ширина диффузионной зоны на два порядка меньше зоны прогрева, поэтому перераспределение легирующих элементов происходит при температуре, равной температуре поверхности, что также удобно для численного исследования, так как позволяет использовать в диффузионной задаче только эту температуру. Исключение составляют ситуации, когда коэффициенты диффузии Б у по порядку сравнимы с коэффициентом температуропроводности к=Х/ср, что для твердой фазы большая редкость. В окрестности температуры плавления это может привести к необходимости построения интерполяционной зависимости температуры от координат с целью использования ее в задаче диф-
фузии, что не представляет большой проблемы. Для представленных в статье расчетных данных в этом не было необходимости.
Обращает на себя внимание качественное и количественное различие в характере распределения концентраций в диффузионной зоне для постоянного коэффициента теплопроводности (штриховая кривая на рис. 5, б, в) и коэффициента X, зависящего от температуры (сплошная кривая на этих же рисунках). Учет перекрестных диффузионных потоков для данного набора параметров оказывается менее существенным (штриховая и штрих-пунктирная кривые).
Нелинейные зависимости концентраций в точке х = 0 от времени (рис. 6, б, в) связаны с нелинейностью задачи даже при постоянной теплопроводности. Так как в модели не учитываются зависимость теплофизических свойств (теплопроводности, плотности, теплоемкости) от концентраций и химические превращения, на поведении температуры появление перекрестных диффузионных потоков не сказывается (см. рис. 5, а и 6, а).
Приведем некоторые цифры, подтверждающие эффективность предложенного алгоритма. При фикси-
Рис. 6. Изменение температуры (а) и концентраций легирующих элементов (б, в) в поверхностном слое в процессе остывания материала с покрытием. Штрих-пунктирная кривая — расчет без учета перекрестных диффузионных потоков и зависимости теплопроводности от температуры, X = 0.45 Вт/(см • с • К); пунктирная кривая — расчет с учетом перекрестных потоков, X = 0.45 Вт/(см • с • К); сплошная кривая — расчет с учетом перекрестных диффузионных потоков и зависимости теплопроводности от температуры. бЦ = 1; = 10-1; Бл = 10-1;
Б02 = 5 • 10 1; Е1 = Е2 = 140 000Дж/моль. Для задачи теплопроводности Ах = 0.001; Аг = 0.00025; для задачи диффузии Дх = 5 • 10-7; Аг = = 1.25 • 10-6
рованных параметрах (указанных в тексте) для расчета одного варианта задачи неизотермической диффузии с плавлением и кристаллизацией на единой разностной сетке (процессор АГНЬ0К-2410) требуется в среднем около часа. При разделении разностных сеток для теплопроводности и диффузии, но при использовании явной схемы, время счета составляет около 10 минут. Последний алгоритм позволил сократить время счета до 1 минуты и снять ограничения (11).
4. О неотрицательности решения задачи диффузии
Известно, что задачи теории многокомпонентной диффузии допускают физически нереальные отрицательные решения, если не выполнены определенные требования. В рамках теории многокомпонентной диффузии коэффициенты диффузии должны удовлетворять соотношениям Онзагера [9], которые в общем случае имеют сложный вид. В случае идеальной смеси удается показать [10], что соотношения Онзагера приводят к выполнению условий, следующих для систем уравнений типа
Ak (u)-
дщ
dt
3 n d
= 22#
i=l l =1 dxi
(
Bkl (u)
du¡
dxi
+ fk (u)
с достаточно гладкими коэффициентами, где и — вектор с компонентами м1, и 2,..., ип. Для нашей трехкомпонентной системы эти условия сводятся к следующему. Для недиагональных коэффициентов должно выполняться Б12 = 0, если а1 = 0; Б21 = 0, если а2 = 0. Следовательно, в (7)
А02 = Б22 = Б52 - 0; б20! = ^ = б24! - 0.
Диагональные элементы матрицы должны быть неотрицательны. Это дает
бЦ + Бла2 + Б^а^ > 0,
ll
lla2
Б22 + Б22а1 + Б22а1 > 0
Эти условия сохраняются и при учете зависимости коэффициентов диффузии от температуры. При постоянстве диффузионных коэффициентов требуется специальный выбор начальных данных [10].
Таким образом, в работе предложена модель формирования диффузионной зоны в процессе остывания материала с покрытием с учетом перекрестных диффузионных потоков, а также алгоритм численного решения подобных задач, учитывающий различие характерных масштабов процессов теплопроводности и диффузии в конденсированных средах. Показано, что использование различных, но согласованных между собой, пространственных сеток для задач теплопроводности и диффузии уменьшает время численного решения в десятки и даже сотни раз по сравнению с традиционными алгоритмами.
При учете образования различных соединений и фаз в процессе диффузии в уравнениях диффузии и теплопроводности появляются дополнительные слагаемые, зависящие от кинетических особенностей физико-химических превращений и приводящие к дополнительным особенностям в численных алгоритмах.
Аналогичный алгоритм удобен и при решении других задач [1, 11].
Работа выполнена при финансовой поддержке РФФИ, грант № 05-08-33412_а.
Литература
1. Букрина Н.В., Князева А.Г. Моделирование неизотермической диф-
фузии в трехкомпонентной системе при остывании наплавленного покрытия // Тезисы докладов 14 Всероссийской школы-конференции молодых ученых «Математическое моделирование в естественных науках», Пермь, Россия, 4-7 октября, 2005 г. - Пермь: Изд-во ПГТУ, 2005. - С. 14.
2. Дураков В.Г., КнязеваА.Г., Прибытков ГА. Формирование диффузионной зоны в подложке в процессе электронно-лучевой наплавки покрытия // Proceedings of 6th International Conference on Modification of Materials with Particle Beams and Plasma Flows. -Томск: Курсив, 2002. - С. 613-616.
3. Тихонов АН., Самарский А.А. Уравнения математической физики. - М.: Наука, 1972. - 735 c.
4. Криштал М.А., Захаров П.Н., Мокров А.П. Диффузия примесей в железоуглеродистых системах // Физика металлов и металловедение. - 1972. - Т. 33. - Вып. 4. - С. 794-799.
5. Карслоу Г., Егер Д. Теплопроводность твердых тел. - М.: Наука,
1964. - 489 c.
6. Князева А.Г., Савицкий А.П. Оценка объемных изменений в диффузионной зоне. 2. Диффузия компонентов в составной пластине конечной толщины // Известия вузов. Физика. - 1997. - Т. 40. -
№ 6. - С. 48-55.
7. Пасконов В.М., Полежаев В.И., Чудов Л.А.Численное моделирование процессов тепло- и массообмена - М:. Наука, 1984. - 288 с.
8. Ворошнин Л.Г., Хусид Б.М. Многокомпонентная диффузия в гетерогенных сплавах. - Минск: Вышейшая школа, 1984. - 142 с.
9. Де Грот С., Мазур П. Неравновесная термодинамика. - М.: Мир,
1964. - 456 с.
10. Вольперт А.И., Посвянский В.С. О положительности решений уравнений многокомпонентной диффузии и химической кинетики // Химическая физика. - 1984. - Т. 3. - № 8. - С. 1200-1205.
11. Букрина Н.В., Князева А.Г., Овчаренко В.Е., Псахье С.Г. Численное исследование формирования переходной зоны между частицами и матрицей в процессе неравновесной электронно-лучевой модификации поверхности композиционного материала // Физ. мезомех. - 2005. - Т. 8. - Спец. выпуск. - С. 53-56.
Приложение
Для задачи многокомпонентной диффузии, описываемой системой n уравнений:
dai
= 2
dt k=1 dx
d ( „ dai ')
D
dx
(A1)
имеем разностную схему вида
„j+1 „ j
At
1
Ax
Dii.i + Dii.i+1
l,i+l
j+1
— a
j+i
l.i
Ax
П1 + П1 а 1+1 а 1+1
ии,1 + Пи,1 -1 л1,1 — л1,г-1
k ф1
Ах
Ах
П1 + П1 а 1+1 — а 1+1
П1к ,г + П1к ,г+1 ае,!'+1 ak,i
2 Ах
П1 + П1 л 1+1 л 1+1
и\к ,i + и\к ,i—1 лl,k лl,k—1
2 Ах
Следовательно, каждое из разностных уравнений диффузии может быть представлено в форме
— Ck,iлk,i + Bk,iлk,i+l = —^^ , (А2)
& = 1, 2, ..., п, г = 1, 2, ..., N где N — число точек разностной сетки.
Система линейных алгебраических уравнений (А2) для каждого диффузионного уравнения может быть без проблем решена независимо от других методом прогонки, причем нетрудно показать справедливость неравенств
|с| > Л + В |, г- = 1, 2, ..., N — 1.
Особенности существуют лишь в аппроксимации граничных условий.
Разложим концентрации а1 в точке х = Ах в ряд Тейлора в окрестности точки х = 0, откуда найдем
(а 2 Л
/ да1 Л дх
а1,1 а1,0
Ах
д а1 дх2
Ах
— + о (Ах 2). (А3)
Следовательно, в окрестности точки х = 0 имеем две системы уравнений:
У *■
к,1 — аk,0
0 k=1
— ( Dik — У Не
Ах Ах I
Т J
( д 2 Л д а*
дх 2 V У 1 0
(А4)
даг-
~д~
V у
=2
0 ¡=1
П
ае,1 ае,0
Ах
Ах ( 2 Л д ае
2 дх2 V У 0 _
(А5)
'дД* Л1 4П/,., — Щ2 — 3П
дх
V J
1*,1 п,*,2
2Ах
; г = 1, 2, ..., п.
Определяя вторые производные концентраций из решения системы линейных уравнений
2
¡=1
(д 2 Л д а*
+ Fi = 0, г = 1, 2, ..., п,
дх 10 0
следующей из (А4), (А5), где
2(уг*к —
ае,1 — ае,0
Ах2
(А6)
и подставляя полученные производные в (А5), найдем разностную аппроксимацию граничных условий.
В частности, для диффузионной задачи (2), (4) имеем
да1 Л а1,1 — а1,0 а2,1 — а2,0
"зТ 1 =^п-------л------+ У12---- -------
дt 10 Ах Ах
— I П11 — Ун -у
— | П12 — У12 “2“
(д 2а Л
а1
дх2
(д2 Л д а2
"дхг
0
да2 Л а1,1 — а1,0 а2,1 — а2,0
ЗТ I =У21 7--------------------- + У22----------------
дt 10 Ах Ах
— 1 П21 —У 21 “2“
(д 2а Л
а1
— 1 П22 — У 22'
Ах
21
(д2 а1 Л дх 2
V У
( д 2 а1 Л
дх 2
V У
+ Л-
( 2 Л д2а
12
дх 2
( 2 Л
дх2
+ Л
22
дх 2
0
+ F2 = 0,
Следовательно,
F2 Л12 F1Л22
д 2 а1 дх 2 0
V /0
(д а2 Л F1Л21 — F2Л11
в
где в = Л22 Л11 Л12 Л21 ■
Арифметические преобразования и разделение переменных, относящихся к разным компонентам, дают
1+1 а1,0
1+1 2,0
(
1 +
А/ Гц Ах к
1+
А/ Г22 Ах к
= ^Гк +1 +
а
Ах к 1
+ п1 + А/ г12 (а ] а] )
+ а1,0 +Т---------Т (а2,1 — а2,0 ) ,
Ах к
= ^Г1 а1+1 +
Ах к 2Д
+ п1 + А/ Г21 (а ] а] )
+ а2,0 + Т----Т (а1,1 — а1,0) ’
Ах к
+
где
Ах
Г11 = П11 — ^ [(Л12 П11 — Л11 П12)(7 21к — П21) + + (Л21П12 — Л22 П11)(Т11к — П11)],
Ах
г12 = П12 —~^ [(Л12 П11 — Л11 П12)(У 22 к — П 22 ) + + (Л21П12 — Л22 П11)(Т 12 к — П12)],
Ах
г21 = П 21 [(Л12 П 21 — Л11П 22)(У 21к — П 21 ) +
+ (Л21 П22 — Л22 П 21 )(711к — П11)],
Ах
г22 = П22 [(Л12П21 — Л11П22)(Т22к — П22) +
+ (Л21П22 — Л22 П21)(^12 к — П12)], что позволяет без труда найти первые прогоночные коэффициенты.
