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

УДК 519.624

И. Ф. Соловьева, доцент ОБ ОСОБЕННОСТЯХ РЕШЕНИЯ ЗАДАЧИ ТРОЕША

In job for the decision of nonlinear regional problems the method of plural bilateral shooting is offered. Practical realisation of a method is considered on an example of the decision of known regional problem Troesh. For the decision of problem Troesh computing schemes of a method of plural bilateral shooting are under construction. They comprise procedure of the decision of problems Cauchy in direct and return directions on shooting subintervals. On an example of problem Troesh it is shown, that the choice of number and lengths of subintervals of shooting provides necessary properties and qualities of problems Cauchy. During the decision of the given problem come to light and properties of matrixes Jakoby for closing systems of the equations are studied. These systems of the equations turn out enough a low order. Borders of a spectrum of matrix Jakoby are defined and their numbers of condi-tionality are calculated. The offered technique allows to solve wide classes of regional problems with an interface, and thus reduces computing difficulties of boundary problems.

Введение. Проблеме численного решения нелинейных двухточечных граничных задач с пограничными слоями уделяется в настоящее время все большее внимание. Эти задачи являются математическими моделями диффузионно-конвективных процессов или родственных физических явлений [1].

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

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

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

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

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

Практическая реализация метода множественной двусторонней пристрелки и его качества зависят главным образом от того, какие

имеются возможности на следующих его основных этапах:

1) выбор точек пристрелки;

2) выбор точек сшива решений;

3) выбор параметров пристрелки и их локализация;

4) выбор длин положительных и отрицательных подынтервалов пристрелки;

5) регулировка свойств замыкающей системы уравнений и ее оптимизация по числу уравнений;

6) организация итерационных процессов и их оптимизация, сведение исходной граничной задачи к совокупности задач Коши. Для решения задач Коши существует много хорошо разработанных методик [2].

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

Математическая модель этой задачи имеет следующий вид:

У1 = У2,

y'2 = kshky1, 0 < t < 1, k > 0 с граничными условиями

уД0) = 0, у,(1) = 1.

(1)

(2)

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

проверки методов решения неустойчивых нелинейных граничных задач. Задача вида (1), (2) получила название задачи Троеша [1].

Решение задачи Троеша. Проведем анализ задачи Троеша.

Матрица Якоби для системы полученных уравнений будет иметь вид

д/ (X, У) ду

0

1

к2 екку1 0

При фиксированном времени матрица Яко-би характеризуется собственными значениями

\2(0) = ±к^екку1(1).

Из последнего выражения следует, что Х1>2 (0) = ±к, \2 (1) = ±^сйк.

Если значения к небольшие, то небольшими будут и значения Л,12(1). Однако, чем больше значение к, тем более сложным является характер поведения решения у1 (X) и его производной у2(^).

Предлагаемый метод множественной двусторонней пристрелки состоит в следующем.

1. Выбираем точки пристрелки:

а = t2<...< ^ 2т 1< Ьт = Ь

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

2. Строим пристрелочные задачи Коши в прямом и обратном направлениях:

и' = /(X, и), X е 4+-1 = ^;-1 < X < ^2;},

и(Х,У2,-1)

X=4 ^ = У21 -1, У2 ¡-1

Я"

V = /(X, V), X е 3--1 = {?21 -1 > X > ?21 -2},

ЧХ, У21 -1)

, = У2/-1, У2 1-1 е , 1 = 1 m,

(3)

(4)

где X2 ■ -1 - точки пристрелки; X 2 ■ - точки сшива решений; У 21-1 - параметры пристрелки.

3. Для полученных пристрелочных задач Коши составляем замыкающую систему уравнений вида

| и(21, У21 -1 ) - ^21, У21+1 ) = 0, 1 = 1 m, к(v(X0, У1), и (X

2т ' У2 т-1 )) = 0.

(5)

4. Переписываем замыкающую систему вида (5) в операторной форме

где

Н (г) = 0,

Н : Ям ^ Ям, N = тп,

г = У , У2 , ..., УТ2т-1)Т.

(6) (7)

Свойства замыкающих систем уравнений (5), (6) зависят от правой части / исходного уравнения, формы граничных условий, области интегрирования [а, Ь], точек пристрелки и^, У 2 -) и траекторий пристрелки и(X, У21-1), v(X, У2;-1). Эти свойства наиболее полно характеризуются матрицами Якоби для соответствующих отображений Н.

Пусть у^) - искомое решение исходной граничной задачи. Введем обозначения

\ * / *Т \ I

У2 , 1 = У(Ь ыХ г = (У , ..., У2т-1 ) .

*Т

(8)

Пусть известны приближения к-го номера к истинным значениям параметров пристрелки

У1(к} =

У1к} У 21}

У3к} =

у1з }

у2?

Тогда будет выполняться замыкающая система вида

и1 (X(0), У1(к)) = 0,

^(2), У1(к)) -и^(2), У3к)) = 0,

V2(X(2), У1(к)) -и2(X(2), У3к)) = 0,

и1 (X(4), У3к)) -1 = 0.

Или в операторной форме она примет следующий вид:

Н (г(к)) = 0,

где

г(к) = (У1(к), у21), У13), У21))Т.

(к К Т

5. Теперь искомое решение у($) исходной граничной задачи представляем формулой

^), у *;-1), X е ./2--1,

y(x) =

и(^ У * у-l), x е 3 2++-l, 1 = 1 т.

(9)

В случае прямого направления матрицы-блоки для пристрелочных задач Коши будут иметь вид

и-.(X) = Ф2 н(и)и2 ;-l(X), = I, Xе У2++-1,

д/( X, и ( X, у21 -1))

Ф 21 -1(и ) =-

и 2 М( X) =

ди

ди ( X, у21 -1) ди

(10)

(11)

Аналогичным образом получаем матрицы-блоки для пристрелочных задач Коши, решаемых в обратном направлении:

V' 1 -1( X) = Ф 21 -1^) У21 -l(X), V1 -1(X21 -1) = I, X е 32--1,

где

Ф 2 -1 =

д/ (Г, УЦ, У 2 ; -,)) ду

V (л ^ У21 -1) ■ 1

У21 -1 С) =-"---, 1 = 1,

(13)

дУ

21 -1

причем

и 22-) = и 21-1^21 ),

V(21-2) = V (г ) 2 1 -1 2 1 -1 2 1 -2

(14)

Будем считать, что для решения замыкающей системы уравнений используется модифицированный метод Ньютона, основанный на аппроксимации матрицы Якоби. Кратко изложим его алгоритм.

1 Шё®} ^=(- Я(г®).

&

(15)

где

2. Построим последовательные приближения

г(к+1) = г(к) + Аг(к), к = 0,1,2,... (16)

3. Найдем поправки А(к) в методе Ньютона:

Аг(к) = (Ьг(к)т, Аг(к)т, ..., Аг«-/, (17)

Н = (И(к), И3к),..., ег, Г = иг(г(к)).

Для замыкающей системы уравнений в задаче Троеша конкретизируем матрицу Якоби и запишем ее в следующем виде:

дН (г(к)) дг

Ф(2к)

)

о(к) о(к)

Введем обозначения Н0 = (дН (г( )))/ дг и образуем матрицу В = НтйН 0.

Вычислим числа обусловленности матрицы Н0. Для этого определим с начала верхнюю и нижнюю грани матрицы В:

(

г(к+1) = В

„(к)

Л

.(к)

к = 0,1, 2,..., г(к) е Я4.

Верхнюю границу Р(В) матрицы В найдем по правилу

Р(В) = Нш г(к) ,к = 0,1,2,... к II

Введем обозначения

С = в(В)Е-В, ю(к) е Я4,

в(С) = Нш ю

к ^ад

(к)

Определим нижнюю границу а(В) матрицы В

а(В) = Р(В) -Р(С).

Итерации векторов г(к) и ^(к) проводились до тех пор, пока не достигалась точность:

г(4+1) - г(4)

ю(*+1) - ю(*)

, 1102

<- 10-5. 2

Искомое значение чисел обусловленности вычислим по формуле

V* Н 0) = .

Р( В)

а(В)'

Методы решения задач Коши, предназначенные для определения блоков и^- и ^(у71), можно видоизменять таким образом, чтобы свойства, характеризующие метод Ньютона, сохранялись.

Заключение. Предложенная методика позволяет решать широкие классы граничных задач с пограничным слоем и при этом уменьшает те вычислительные трудности, которые присущи традиционным методикам. Система уравнений вида (5) не будет изменяться также в случае неравномерности выбранной сетки. На систему не окажет особого влияния и перемена методов решения задач Коши. А для решения задач Коши в настоящее время существует достаточно большой арсенал хорошо разработанных методик [2].

Литература

1. Холл, Д. Современные численные методы решения обыкновенных дифференциальных уравнений / Д. Холл, Д. Уатт; пер. с англ. - М.: Мир, 1979. - 312 с.

2. Деккер, К. Устойчивость методов Рунге-Кутты для жестких нелинейных дифференциальных уравнений / К. Деккер, Я. Вервер; пер. с англ. - М.: Мир, 1983. - 200 с.

МЕХАНИКА

УДК 531.19

Р. Н. Ласовский, ассистент; Г. С. Бокун, доцент; В. С. Вихренко, профессор КИНЕТИКА ФАЗОВЫХ РАССЛОЕНИЙ В РЕШЕТОЧНЫХ ФЛЮИДАХ

The equations which describe phase layering kinetics in the 2D and 3D lattice fluids are obtained. The evolution of the system with step-like initial concentration distribution is investigated. It is shown that the second condense area appears in the 2D system during evolution. This area disappears with time. The inoculating area critical size is determined. Several complementary condense areas origin in 3D system is observed. Number of areas depends on initial conditions. However, as opposed to 2D system, phase layering in 3D system takes place with any nucleator size.

Введение. Описание неоднородных неравновесных состояний систем многих частиц и сопутствующих им процессов является одной из наиболее сложных и важных проблем статистической механики. Процесс зарождения новой фазы, структурные превращения в системе относятся к числу вышеотмеченных проблем, и их рассмотрению посвящено много работ [1-4], использующих, как правило, классическую теорию нуклеации или прямое компьютерное моделирование и относящихся, в основном, к жидкостям и газам.

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

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

Эта модель используется в тех случаях, когда рассматриваемая система может быть представлена двумя подсистемами: относительно жесткой несущей подсистемой, создающей некоторый потенциальный рельеф, и подсистемой лабильных частиц. Примерами таких систем являются водород в металлах [8], литий в оксидах металлов [9] или в углеродных матрицах [10].

1. Дифференциально-разностное уравнение кинетики межфазной границы. Как показано в работе [11], эволюцию плотности частиц

интеркалянта можно описать уравнением неразрывности

^ = j < I dt j »

>.

(i)

где суммирование ведется по всем первым соседям узла 7.

Средний поток через границу между узлами 7 и ] в уравнении (1) определяется следующим выражением:

< 1Ц > = ^ ^(0г ,0; )(ер^ - е^), (2)

где ^о - частота перескоков частицы в пределе низких концентраций; ^(07, 0;) - неравновесная функция распределения, определяющая вероятность того, что ближайшие соседние узлы 7 и ] вакантны; р = (квТ)-1 - обратная температура;

кв - постоянная Больцмана; ц7 - локальное значение химического потенциала, отнесенное к узлу 7. Бинарная функция распределения имеет вид

(3)

где р07 - концентрация вакансий в 7-м узле.

Химический потенциал определим выражением

= Pi»

P

Пп,

(4)

0i jФ»

где р17 и р07 - концентрация частиц или вакансий в 7-м узле соответственно; 'Цji - величина, рассчитываемая как положительный корень уравнения

п,, +п

Pi, -Po, - W(Pi, -Pi,) Wpi

= 0, (5)

Po, Po,

в котором W = exp (-pj), J - энергия взаимодействия ближайших соседей.

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

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

При решении системы кинетических уравне-иий приняты периодические граничные условия.

0,7

§ 0,6 и

^ 0,5

о Я

ё 0,4 У

0,3 0,2

Рис. 1. Зарождение конденсированной фазы в двухмерной системе. Время указано в шагах алгоритма Эйлера: 1 - 0; 2 - 5 -103; 3 - 5-105;

4 - 2,5-106; 5 - 5-106

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

Кроме того, на временах, соответствующих «104 шагам алгоритма Эйлера, зарождается

вторая область на границе ячейки, в которой до « 2,5 -105 шагов концентрация постепенно увеличивается. Это явление обусловлено конкуренцией межчастичного притяжения, термоак-тивационного механизма перемещения частиц и наличия вакансий. Однако центральная часть с течением времени продолжает доминировать, и система расслаивается на две четко выраженные фазы, разделенные переходными слоями, концентрации в которых соответствуют равновесным.

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

На рис. 2 показана кинетика распределения концентрационного поля в трехмерной системе, в которой начальное распределение также имеет ступенчатую «затравку» в центре. Неоднородность концентрации задавалась в двух измерениях и в вертикальном направлении удерживалась постоянной. В горизонтальной плоскости использовались периодические граничные условия.

Через « 6 • 104 шагов в системе начинают появляться восемь конденсированных областей на гранях системы и одна диагональная. Концентрация в этих областях растет до «9 -104 шагов, потом к «2 -105 шагам области на гранях, притягиваясь попарно друг к другу, сливаются и образуют четыре пика, которые исчезают к «2,5 -105 шагам. Механизм этого процесса схож с вышеописанным механизмом, происходящим в двухмерной системе. Таким образом, в системе остаются две конденсированные фазы: центральная и диагональная.

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

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

1 ' Г-1— 1 1 -»-1 ■ 1

- И 3 и -

\г4

- \ / 2\ / -

1

1 ■ |_|_ ■ 1 _|_1 ■ '

0 20 40 60 80 100 Номер слоя