2.6. Методы решения задач с несколькими пространственными переменными (Математическое моделирование задач газодинамики и тепло-массообмена. Молчанов А.М. 2013)

6. Методырешениязадачснесколькимипространственными переменными6.1.Особенности двух- и трехмерных задач.В задачах с несколькими пространственными переменными последискретизации получается гораздо более сложная система алгебраическихуравнений, чем для задачи с одной пространственной переменной.Рассмотрим, например, частный случай обобщенного уравнения (1.16),когда задача является двумерной, а конвекцией и источниковым членомможно пренебречь. Полагаем для простоты, что Γ = const , ρ = const .

В этомслучае уравнение имеет вид∂Φ∂ 2Φ∂ 2Φ= ν 2 +ν 2∂t∂x∂yЭто уравнение(6.1)можно интерпретировать и как дифференциальноеуравнение теплопроводности. В этом случае под Φ подразумеваетсятемпература T , а под ν - коэффициент температуропроводности.Пусть задано начальное распределение Φ и граничные условия первогорода:Φ ( 0, x, y ) = φ0 ( x, y ) ,Φ w = ψ ( t , x, y )(6.2)Построение разностных схем для многомерных задач аналогично этомупроцессу для одномерной задачи. Получение многомерных дискретныханалогов проводится с помощью прямого распространения результатов,полученных при рассмотрении одномерной задачи.Используем метод, основанный на разложении в ряд Тейлора. Неявныйразностный аналог задачи (6.1) на сетке, изображенной на рис.5, имеет вид:Φ nm+, k1 − Φ mn , k−νΦ mn ++11, k − 2Φ nm+, k1 + Φ nm+−11, k∆t∆x 2n = 0,1, 2,3,..., Nτ ; Nτ ∆t = tmaxm = 1, 2,3,..., N X ;k = 1, 2,3,..., NY ;−νΦ nm+, k1 +1 − 2Φ nm+, k1 + Φ nm+, k1 −1∆y 2= 0;,(6.3)( N X − 1) ∆x = LX ;( NY − 1) ∆x = LYгде LX , LY - размеры расчетной области по осям x, y ,Φ nm , j - решение разностной задачив момент времени t n = n∆t , в точке( m, k ) , которая находится на пересечении линий(см.

рис. 5).Рис.5. Сетка для двумерной задачи.xm = ( m − 1) ∆x и yk = ( k − 1) ∆yУравнение (6.3) приводится к видуb Φ nm+, k1 +1 + a Φ nm+, k1 + c Φ mn +, k1 −1 + d Φ mn ++11, k + e Φ mn +−11, k = f ,(6.4)гдеa =1+ 2ν∆t+2ν∆t,∆x∆y 2ν∆tν∆tb=c=− 2, d =e=− 2,∆y∆xУравнение (6.4)2(6.5)f =Φnm, jсодержит 5 неизвестных. Вся система линейныхалгебраических уравнений описывается так называемой пятидиагональнойматрицей коэффициентов.Аналогично для трехмерной задачи получается семидиагональнаяматрица.В отличие от одномерной задачи, которая описывается трехдиагональнойматрицей коэффициентов и которая эффективнорешаетсяметодомпрогонки, для задач с несколькими пространственными переменными такихэффективных методов нет.Прямые методы решения алгебраических уравнений,применяемые кдву- или трехмерным задачам, становятся более сложными и требуютсущественно большей машинной памяти и затрат вычислительного времени.Альтернативойявляютсяитерационныеметодыдлярешенияалгебраических уравнений, а также подход построения разностных схем,основанный на методах расщепления.УпражнениеОпределить порядок аппроксимации разностной схемы (6.3).УпражнениеДоказать, что для разностной схемы (6.3) спектральный признак Нейманавыполняется при любых значениях шагов по времени и по пространству.Указание.

Для решения этой задачи необходимо доказать ограниченностирешений вида Φ nm, k = λ nei (α m + β k ) при любых значениях шагов.6.2.ИтерационныйметодГаусса-ЗейделяспеременойнаправленийПростейшим из всех итерационных методов является метод ГауссаЗейделя, в котором значения переменной рассчитываются путем обращения вопределенном порядке к каждой узловой точке. В памяти вычислительноймашины держится только один массив значений Φ m, k . По мере обращения кочередной узловой точке соответствующее значениеΦ m, kв памятивычислительной машины (начальное приближение или значение Φ m, k спредыдущей итерации) заменяется на новое.Для решения системы (6.4) на каждом шаге по времени используются двапрохода в направлении основного движения потока (в данном случае x )Для j = 1,3,...

( j - номер итерации):1) Проход в обратном направлении ( m = N X − 1, N X − 2,...,3, 2 )b Φ (mj,)k +1 + a Φ (mj,)k + c Φ (mj,)k −1 + d Φ (mj+) 1, k + e Φ (mj−−1,1)k = f(6.6)Для m = N X − 1 значение Φ (mj+) 1, k = Φ (Nxj ), k определяется граничными условиямина правой границе и, таким образом, является известной величиной. Врезультате при m = N X − 1 получается система, которая в каждой строкесодержит только 3 неизвестные величины - Φ (mj,)k −1 , Φ (mj,)k , Φ (mj,)k +1 . Таким образом,система описывается трехдиагональнойэффективно решается методом прогонки.матрицейкоэффициентовиДалее, при m = N X − 2 используется уже найденное значение Φ (mj+) 1, k = Φ (Nxj )−1, k ,и мы снова получаем трехдиагональную матрицу.Этот процесс продолжается до слоя m = 2 , а далее начинается проход вдругом направлении.2) Проход в прямом направлении m = 2,3,...N X − 2 , N X −1b Φ (mj,+k1+)1 + a Φ (mj,+k1) + c Φ (mj,+k1−)1 + d Φ (mj+) 1, k + e Φ (mj−+1,1)k = f(6.7)Для m = 2 значение Φ (mj−+1,1)k = Φ1,( jk+1) определяется граничными условиями налевой границе и является известной величиной.

Решение полученнойсистемы, а также систем при всех последующих значений m , полностьюаналогично предыдущему пункту.Отметим, что для простоты записи в формулах (6.6) и (6.7) приобозначении Φ не указывается номер шага по времени – везде используется( n + 1) -ый номер.Для начальной итерации используется значение с предыдущего шага повремени:Φ (m,)k = Φ nm , k0(6.8)Для повышения эффективности метода можно чередовать направлениепроходов, а также использовать проходы не только в продольномнаправлении, но в поперечном.Главным недостатком метода Гаусса-Зейделя является медленнаясходимость, особенно когда используется большое число узловых точек.6.3.Метод приближенной факторизацииМетоды факторизации состоят в том, что вместо решения одной системылинейных алгебраических уравнений решаются несколько систем, каждая изкоторых позволяет использовать более простые методы решения.Систему (6.4) можно записать в обычном векторном видеAx = b ,(6.9)a c .

. e . . . . . . b a c . . e . . . . . . b a c . . e . . . .. . b a c . . e . . .d . . b a c . . e . . A= . d . . b a c . . e ., . . d . . b a c . . e. . . d . . b a c . .. . . . d . . b a c . . . . . . d . . b a c . . . . . . d . . b a(6.10)где. n +1  Φ m +1, k .. n +1  Φ m , k +1 x =  Φ nm+, k1  , Φ nm+, k1 −1 .. Φ nm+−11, k .. nΦ m +1, k .. n  Φ m , k +1 b = Φ nm , k Φ nm , k −1 ..Φ nm −1, k .(6.11)Матрица A имеет N строк и N столбцов, где N = N X × NY . Соответственноразмерность векторов x и b также равна N .Метод приближенной факторизации состоит в том, что матрицапредставляется в видеA ≅ AX ⋅ AY ,(6.12)где ax .dAX =  . . . .Здесь.ax.e.ax.e...e...d.....d..ax.d..ax.de.ax.коэффициенты ayb ..

 , AY =  . .e. . .a x ...b,c ,d ,ecayb.....cayb...совпадают..cayb..с...cayb.....cayb. . .. . ca y (6.13)соответствующимикоэффициентами матрицы A , а коэффициенты a x и a y зависит только отшагов ∆x и ∆y соответственно:ax = 1 + 2ν∆t∆x 2ν∆tay = 1 + 2 2∆y(6.14)Вводим обозначениеx = AY x(6.15)AX x = b(6.16)и решаем системуЕсли подробно расписать это выражение, то получим систему в виде m ,k + d Φ m +1,k + e Φ m −1, k = fax Φ(6.17)Эта система имеет трехдиагональную матрицу коэффициентов иэффективно решается методом прогонки.Затем решается система (6.15) m,k ,b Φ mn +,1k +1 + a y Φ nm+,1k + c Φ nm+,1k −1 = Φкоторая тоже имеет трехдиагональную матрицу коэффициентов.(6.18)Основным недостатком метода приближенной факторизации является то,что ни один из ненулевых элементов матрицыAне совпадает ссоответствующими элементами матрицы AX ⋅ AY , а некоторые нулевые дажестановятся ненулевыми.

Это может приводить к большой погрешностиметода и медленной сходимости.Для трехмерных задач метод приближенной факторизации сходитсякрайне медленно и его не рекомендуется использовать.6.4.Модифицированный метод приближенной факторизацииМаккормака [4]Метод приближенной факторизации может быть улучшен с тем, чтобыуменьшить ошибку разложения матриц. Рассмотренный далее методпозволяетполучитьточноесовпадениененулевыхэлементовфакторизованной и исходной матриц.

Правда, некоторые нулевые элементыисходнойматрицыстановятсяненулевымиприфакторизации,новозникающая ошибка может быть ликвидирована с помощью итераций.Матрица A , входящая в систему (6.9), представляется в видеA ≅ AX D −1 AY ,(6.19).a c . . . . . . .b a c . . . .  . b a c . . .. .e .  , AY =  . . b a c . .  . . . b a c .. ea . .

. . . b a c . . . . . b a. a (6.20)гдеa.dAX =  .....e..a.e..a.ed..a..d.a..d....d.Здесь коэффициенты все a, b , c , d , eсовпадают с соответствующимикоэффициентами матрицы A .Диагональная матрица D состоит только из элементов a :a . . . . . . . a . . . . .. . a . . . .D= . . . a .

. .. . . . a . .. . . . . a . . . . . . . a(6.21)Таким образом, вместо системы (6.9) решается системаAX D −1 AY x = b(6.22)Решение этой системы не представляет большой сложности: оно состоитиз обращений двух трехдиагональных матриц и одного умножения надиагональную матрицу.Для уменьшения погрешности этого алгоритма можно добавитьследующий итеративный процесс.Систему (6.22) можно представить в виде:AX D −1 AY x = Ax + Px ≡ b + Px ,(6.23)где Px - невязка – отклонение приближенного решения от точного,которое находится из решения системыAx = bМожно ввести итеративный процесс решения основнойосновываясь на уравнении (6.23):AX D −1 AY x( j ) = b + Px ( j −1) ,j = 1, 2,....где j - номер итерации.Для нулевой итерации:(6.24)системы,(6.25)x ( ) = [ 0]0(6.26)После первой итерации, которая фактически описана выше (формула(6.22) ), необходимо найти невязку.

Из формулы (6.23) следует:Px (1) = AX D −1 AY x(1) − Ax (1) = b − Ax(1)(6.27)Затем решаем систему (6.25) для второй итерации j = 2 . И т.д.Для последующих итераций имеем:Px ( j ) = AX D −1 AY x ( j ) − Ax( j ) = b + Px( j −1) − Ax( j ) ,(6.28)Если ввести обозначение для правой части (6.25)r ( j −1) = b + Px( j −1)(6.29)то невязка после j -ого шага определяется по формуле:Px ( j ) = r ( j −1) − Ax( j )(6.30)Для хорошей сходимости достаточно двух итераций.В работе [4] показано, что метод модифицированной факторизациинамного эффективнее метода приближенной факторизации как для решениядвумерных задач, так и для трехмерных.6.5.Методы расщепленияМетоды, описанные в предыдущих разделах, основаны на итерационномрешении системы (6.4), которая была получена при использованииразностной схемы (6.3) для задачи (6.1).Разностные схемы расщепления основаны на ином подходе – в них самиразностные схемы строятся так, что их решение не представляет большойсложности.Основную идею схем расщепления можно изложить так.Рассмотрим дифференциальную задачу вида∂Φ= AΦ, ∂tΦ t = 0 = φ0 (6.31)где A - некоторый дифференциальный оператор по пространственнымпеременным.Например, для задачи (6.1) этоAΦ = ν∂ 2Φ∂ 2Φ+ν∂x 2∂y 2(6.32)Допустим, что правая часть уравнения (6.31) имеет видAΦ = A1Φ + A 2 Φ(6.33)На каждом шаге по времени при переходе от момента времени t = t n кt = t n +1 = t n + ∆t заменим решение задачи (6.31) на решение двух задач:∂v= A1v, t n ≤ t ≤ t n +1 , ∂tv ( t n , x, y ) = Φ ( t n , x, y ) (6.34)∂w= A 2 w, t n ≤ t ≤ t n +1 ,∂tnn +1w ( t , x, y ) = v ( t , x, y ) (6.35)w ( t n +1 , x, y ) = Φ ( t n +1 , x, y ) + 0 ( ∆t 2 ) ,(6.36)иМожно показать [2], чтот.е.