9.1.3, Схема чехарда и схема Лакса - Вендроффа

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

2а/ + ja. -О (9.14)

и соответствующему ей алгоритму

Формула (9.15) соответствует уравнению (9.2) с ошибкой аппроксимации, имеющей порядок 0(Д/2, Дл:). Применение к (9.15) анализа устойчивости по Нейману дает коэффициент усиления G, представленный в табл. 9.1. Как можно видеть, G = 1.0, если С 1, однако при некоторых значениях 6 получим G > 1.0, если С > 1. Следовательно, если удовлетворяется условие КФЛ, схема чехарда будет нейтрально устойчивой. Это и желательно, так как в этом случае решение не затухает с течением времени; отметим, что на то же самое указывает и физический характер процесса конвекции согласно (9.2).

Как нетрудно видеть, решение обеспечиваемое фор-

мулой (9.15), как бы перепрыгивает через решение Г ; отсюда и название -схема чехарда . Однако структура формулы (9.15) указывает на то, что вычислительная сетка уподобляется здесь шахматной доске; решение, строящееся для черных полей , оказывается совершенно независимым от решения для белых полей . Это обстоятельство может приводить к появлению двух расщепившихся решений. Роуч [Roache, 1972] обсуждает различные варианты стратегии для преодоления данного затруднения, одним из которых является осреднение решения по двум временным слоям, а именно п и n-f 1.

Халтнер и Уилльямс [Haltiner, Williams, 1980] провели анализ применения схемы чехарда (9.15) к задаче с гармоническим начальным условием, т. е. когда функция в правой части формулы (9.4) принимает форму f (х)= ехр(г/тДл:), где т - волновое число (§ 9.2). Они получили для (9.15) аналитическое решение, состоящее из двух волн или двух мод, происхождение которых связано с тем, что формула (9.15) коррелирует между собой три временных слоя. Одна из этих мод соответствует физическому решению, а вторая - это паразитная вычислительная мода, распространяющаяся в направлении.

противоположном моде физического решения, причем величина этой моды меняет свой знак на каждом шаге по времени.

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

Схема чехарда является трехслойной, и поэтому при ее применении требуется некий альтернативный подход для реализации первого шага по времени. Эту схему рекомендуется применять [Morton, 1971] для решения гиперболических задач, например задач, описываемых уравнением конвекции (9.2); заметим, что эквивалетная схема Ричардсона для уравнения диффузии является неустойчивой.

Схема Лакса - Вендроффа представляет собой весьма эффективный (и популярный) алгоритм для решения уравнений, описывающих течение невязкой сжимаемой жидкости (§ 10.2 и п. 11.6.1). В применении к уравнению конвекции (9.2) схема Лакса - Вендроффа совпадает с методом Лейза [Leith] (см. [Roache, 1972]). Оба метода основаны на построении приближения второго порядка для члена с помощью разложения в ряд Тейлора, в результате которого получается представление

После введения представления для дТ/дх с применением центральных разностей мы приходим к конечно-разностной модели для уравнения (9.2) по Лаксу - Вендроффу в форме тп Q 5 j,n Q + rf+i). (9.16)

Схема Лакса - Вендроффа согласуется с уравнением (9.2) при ошибке аппроксимации порядка 0(Л2, Ах) и является устойчивой, если удовлетворяется условие КФЛ С 1.0 (табл. 9.1). Численные решения с использованием схемы Лакса - Вендроффа приводятся в п. 9.1.5.

9.1.4. Схемы Кранка -Николсона

Схема Кранка - Николсона оказывается весьма эффективной в применении к одномерному уравнению диффузии. В данном разделе конечно-разностные и конечно-элементные схемы

Idt if

в которую входит массовый оператор Л! = {-,-f, }. Введение разностного представления dT/dt по Кранку - Николсону приводит к схеме

Мх + uLx [0.5 {Tf + Т-г)] - О, (9.20)

которую можно сравнить со схемой (9.17). Массовый оператор можно представить в более общем виде, аналогичном формуле (8.44), а именно

Л1,{б, 1-26, 6}, (9.21)

где значение 6=1/6 дает конечно-элементную форму, а значение 6 = 0 - конечно-разностную форму. После подстановки (9.21) в уравнение (9.20) получается трехдиагональный алгоритм

(6 ~ 0.25С) да -f (1 - 26) -f (6 -Ь 0.25С) Tjtl =

= (6 -f 0.25С)I -f (1 - 26) Ti -f (б - 0.25C) (9.22)

Реализацию этого алгоритма дает программа TRAN (рис. 9.8).

Конечно-элементная форма (9.20) обладает ошибкой аппроксимации порядка 0(А2, Ах). Следовательно, при достаточно

Кранка - Николсона применяются к одномерному уравнению конвекции (9.2).

Конечно-разностная схема Кранка - Николсона может быть записана в виде

7 + (0.5L.r + О.ЬитГ) = О, (9.17)

где LxTj={Tj+i - Ti-i)/2Ax. Из уравнения (9.17) получается следующий трехдиагональный алгоритм:

-0.25С7 + + 7 + + 0.25СГ +/ = 0.25СГ -1 + Т - 0.2507 +

(9.18)

для которого можно построить эффективное решение с помощью алгоритма Томаса (п. 6.2.2). Схема Кранка - Николсона согласуется с уравнением (9.2) и обладает погрешностью порядка 0(А2, Дл:2). Анализ устойчивости по Нейману дает коэффициент усиления G, показанный в табл. 9.1. Ясно, что схема Кранка - Николсона является безусловно устойчивой.

Если к уравнению (9.2) применить метод конечных элементов Галёркина с нелинейной интерполяцией по пространству, то получается следующая система обыкновенных дифференциальных уравнений:

].-f а1,Г; = 0, (9.19)