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

Вычислительные технологии

Том 8, № 3, 2003

ОБ ОДНОМ МЕТОДЕ ПОСТРОЕНИЯ СХЕМ ТОЧНОЙ ФАКТОРИЗАЦИИ ДЛЯ ЧИСЛЕННОГО РЕШЕНИЯ ГИПЕРБОЛИЧЕСКИХ УРАВНЕНИЙ*

Т. З. ИсмАгилов Новосибирский государственный университет, Россия

В. М. КОВЕНЯ

Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: [email protected]

The method for obtaining all finite difference schemes with precise factorization for numerical solution of hyperbolic equations is suggested. The existence of the multi parametrical family of exact factorization schemes is proved. This family schemes can be easily checked for scalar resolvability. The exact factorization schemes for the gas dynamics equations in different various gas dynamics variables are considered as the application.

Введение

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

К настоящему времени для численного решения гиперболических уравнений разработано большое число явных и неявных разностных схем [1-12]. Неявные схемы обычно основаны на методе приближенной факторизации или методе расщепления, позволяющих свести решение исходных задач к последовательности решения их одномерных аналогов или более простых задач. Однако решение многомерных задач связано с определенными проблемами, так как реализация явных схем приводит к ограничениям на устойчивость и, как следствие, к большому числу арифметических операций, а для неявных схем решение сводится к векторным прогонкам, требующим обращения матриц размера m х m в каждом узле расчетной сетки, где m — число уравнений. С увеличением числа уравнений число арифметических операций на узел сетки возрастает по степенному закону как ~ тъ . Для исключения векторных прогонок при решении уравнений газовой динамики в работах [5, 6] предложено расщепление уравнений по физическим процессам, что позволило

* Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований, грант № 02-01-01029, и Минобразования РФ, грант Е 02-1.0-25.

© Т. З. Исмагилов, В.М. Ковеня, 2003.

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

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

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

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

1. Факторизованные схемы

Рассмотрим систему гиперболических уравнений

д£ д£

ж+ВдХ = 0 (11)

где £ — вектор размерности т, а матрица В имеет действительные собственные числа. д

Пусть Л = ——+ 0(Нк) — разностный оператор, аппроксимирующий первую производную дх

с порядком к. Разностная схема с весами

£«,+ 1 _ £п

+ ВпЛ [а£п+1 + (1 _ а)£п] = О

или эквивалентная ей схема в каноническом виде

fn+1 _ f n

(I + TaBn Л)-= _ВПЛГ (1.2)

т

аппроксимируют исходные уравнения с порядком 0(тт + ), где m =2 при а = 0.5 + 0(т). Ее решение может быть получено векторной прогонкой, требующей обращения матриц размера m х m. Для B = const схема (1.2) безусловно устойчива при а > 0.5. Пусть

B = Bi + B2. (1.3)

Приближенно факторизуя оператор I + TaB"^ « (I + тaBnЛ) (I + TaBnЛ) , рассмотрим разностную схему

f n+1 _ fn

(I + TaB^) (I + тaBnЛ)-= _B^fn, (1.4)

т

аппроксимирующую (1.1) с тем же порядком, что и базовая схема (1.2). Очевидно, разностная схема

= _B^fn,

(I + та^Л) Г+1/2 = Г, (1.5)

(I + TaB2^) |n+1 = |n+1/2,

f n+1 = fn + т ^n+1

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

П (/ + таВпА) = I + таВпЛ + т2а2В?АВ^А ¿=1

дополнительных членов Q = т2а2В1ЛВ2Л, т.е. членов второго порядка малости, играющих роль диссипативных членов и отсутствующих в исходной нефакторизованной схеме (1.2). Об их негативной роли говорилось выше. Сформулируем следующую задачу

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

В + В2 = в, В1В2 = 0. (1.6)

Случай, когда В1 или В2 совпадает с В или нулевой матрицей, как тождественный не рассматривается. Тогда устойчивость и свойства факторизованной схемы совпадают со свойствами базовой схемы (1.2).

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

2

2. Схемы с точной факторизацией

Из определения гиперболичности системы (1.1) следует, что матрица В может быть представлена в виде

( ¿1 О О \

О ¿2 О

в = лвь, в

(2.1)

V О О ¿т /

Здесь В — диагональная матрица с различными коэффициентами dj•, а матрица Л = Ь 1, т.е. обратная к Ь. Тогда

В1 = ЬВ1Л, В2 = ЬВ2Л и, очевидно, что условия (1.6) эквивалентны условиям

В + В2 = В, В1В2 = О, (2.2)

т.е. определители В1 и В2 должны быть равны нулю. Таким образом, решив задачу (2.2) для диагональной матрицы В размера т с произвольными и различными коэффициентами, мы решим задачу о нахождении точных факторизаций для любой одномерной гиперболической системы размерности т. После нахождения матриц Dj для ] = 1, 2 значения матриц Bj определяются по формулам

В1 = В2 = ЛВ2Ь. (2.3)

В1В = В2. (2.4)

Опишем алгоритм нахождения матриц Dj•. Исключая В2 из второго уравнения (2.2), получим

В1

Умножая уравнение (2.4) справа на единичный вектор ej, получим (с учетом перестановочности В и В1)

¿1(В1 е, ) = В^е,). (2.5)

Это означает, что или вектор D1ej является собственным вектором В1 или В1 равен нулю, т.е. ]-й столбец матрицы В1 состоит из одних нулей.

Рассмотрим случай, когда размерность системы равна трем. Пусть все столбцы матрицы В1 ненулевые. Тогда у В1 существуют три различных собственных значения ¿г (г = 1, 2, 3). Это противоречит тому факту, что определитель матрицы В1 равен нулю, значит, у В1 есть хотя бы один нулевой столбец. Допустим, что у В1 два нулевых столбца, т.е. существует одно ненулевое собственное значение. Тогда, решая уравнение (2.4) с помощью метода неопределенных коэффициентов, для В1 получаем

/ ОО \ / О а О \ / ОО а \

= а О О , В2 = О ¿2 О , И О О в , (2.6)

\ в о о у \ О в о ) \О О ¿3 /

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

/ О а в \

Вт = | О ¿11 ¿12 0 ¿21 ¿22

Подставляя эту матрицу в уравнение (2.4), получаем, что для матрицы Т должно выполняться равенство

' ¿2 0 о 4

ТТ = Т

(2.7)

где Т = (¿¿-). Определитель матрицы Т равен произведению двух ненулевых собственных чисел и поэтому у Т есть обратная матрица. Умножая (2.7) на матрицу, обратную к Т, слева получаем

' ¿2 0 о 4

Т

т.е. матрица равна

А

0 а в 0 4 0 0 0 4

Повторяя эти рассуждения для других случаев (нулевой второй или третий столбцы), получим оставшиеся решения уравнения (2.4), т.е. значения матриц Я (^ = 4, 5, 6) равны

¿1 0 0 0 ¿2 0

а в 0

¿1 0 0 а0в 0 0 ¿3

Я?

0 а в 0 ¿2 0 0 0 ¿3

(2.8)

Таким образом, пространство решений уравнений (2.4) состоит из шести двухпарамет-рических семейств, которые могут быть объединены в две группы. Первая группа состоит из матриц, в которых один столбец ненулевой (2.6). Вторая группа состоит из матриц, в которых два ненулевых столбца (2.8).

Для системы (2.1) размерности, равной 4, будут 14 семейств решений, которые могут быть объединены в три группы. Первая группа состоит из матриц, у которых один столбец нулевой. В этой группе четыре трехпараметрических семейства матриц Я, причем ]-е семейство получается из матрицы Я заменой ^ на 0 и заменой трех нулей в ]-й строке на три произвольных параметра. Например, для

( 0 а в я \

0 ¿2 0 0

0 0 ¿3 0

\ 0 0 0 )

(2.9)

Оставшиеся матрицы Я (/ = 2, 3, 4) получаются заменой ¿г на нуль и нулей в /-й строке на три произвольных параметра:

я2

/ ¿1 0 0 0 \

а 0 в 7

0 0 ¿3 0

\ 0 0 0 ¿4 )

Я3

(¿1 0 0 0

0 ¿2 0 0

а в 0 7

0 0 0 ¿4 )

я4

( ¿1 0 0 0 \

0 ¿2 0 0

0 0 ¿3 0

V а в 7 0 у

(2.10)

Во второй группе будут шесть параметрических семейств матриц Я, у которых два столбца нулевые. Для построения этого семейства необходимо выбрать на ненулевых столбцах на диагонали элементы ^ и ф, а на пересечении этих столбцов с оставшимися двумя строками поместить четыре произвольных элемента (параметра). Например, для первых двух

нулевых столбцов

D5

0 0 а в ^

0 0 5 Y

0 0 d3 0

0 0 0 d4 )

Оставшиеся матрицы (/ = 6,..., 10) строятся аналогично.

Третья группа состоит из матриц, три столбца которых нулевые. В этой группе семейство из четырех трехпараметрических матриц . Каждое ]-е семейство получается из нулевой матрицы заменой нулевого элемента, расположенного на пересечении ]-го столбца и ]-й строки на dj и заменой трех нулей в ^'-м столбце на три произвольных параметра. Например, для первого ненулевого столбца эта матрица имеет вид

D11

( di 0 0 0 \

а 0 0 0

в 0 0 0

\ 5 0 0 0}

а остальные для I = 11,..., 14 строятся аналогично. Изложенный выше алгоритм может быть обобщен и для системы уравнений более высоких порядков: пусть дана система гиперболических уравнений размерности т. Невырожденное преобразование (2.1) приводит ее к диагональной форме. Выполнение условий точной факторизации (2.2) приводит к системе уравнений (2.4), решение которой и даст все точные факторизации оператора В по формулам (2.3). Тогда факторизованная разностная схема (1.4) или (1.5) совпадает с нефакторизованной схемой (1.2), но ее реализация может оказаться проще для матриц Bj, а выбором параметров можно управлять при построении разностных схем, реализуемых скалярными прогонками.

3. Скалярная разрешимость разностных схем

Пусть имеется система разностных уравнений

(I + таСПЛ) | = {,

(3.1)

где

6

Cm

C11 . .. C1 m \

C = C21 . .. C2m

Cm1 . .. Cmm /

fi

fm

получены из нефакторизованной схемы (1.2) или схемы (1.5) на дробных шагах. Будем считать, что коэффициенты Cj заданы на n-м слое, т.е. известны, и тогда система уравнений (3.1) линейна. Для уравнений с постоянными коэффициентами при Cj = const система уравнений (3.1) в силу гиперболичности исходных уравнений (1.1) может быть приведена к диагональной форме

(I + таВЛ) n = g, (3.2)

где D — диагональная матрица с элементами dj, j = 1, ...,m. Действительно, так как

C = R-1 DR, I = R-1R,

уравнения (3.1) представляются в виде

R-1 (I + та£Л) R£ = f,

эквивалентном (3.2) при n = R£, g = Rf• Очевидно, система уравнений (3.2) может быть решена независимо для каждой компоненты n скалярными прогонками или по схеме бегущего счета в зависимости от оператора Л.

Для уравнений с переменными коэффициентами такое представление невозможно, кроме случая, когда C — диагональная матрица. Предположим, что за счет выбора вектора искомых функций в исходных уравнениях (1.1) или расщепления матрицы B на Bj матрица Cп имеет такую структуру, что некоторые из ее диагональных элементов равны нулю, т.е. 0ц = 0. Тогда

m— 1

& = f - та^сЖг, (3.3)

i=i

где i = /. Если путем подстановок компонент & вектора ^ в другие уравнения из (3.1) система уравнений (3.1) может быть приведена к подсистеме (подсистемам) уравнений, разрешенных относительно отдельных компонент &, т.е. к виду

П& = dl, (3.4)

где П — разностный оператор, полученный в уравнении для & после исключения других компонент, то такую систему уравнений назовем скалярно разрешимой. Поясним вышесказанное на примере. Пусть (3.1) — система уравнений размерности m = 2. Очевидно, среди всех представлений матрицы C существуют четыре формы матриц

Ci = ( ci1 0 ) , C2 = ( ci1 ci2 ) , Сз = ( 0 ci2 ) , C4 = ( ci1 ci2 ) , V C12 C21 J V 0 C12 / V C2i C22 ) V C2i 0 У

при которых разностная схема реализуется скалярными прогонками. Конечно, отдельные компоненты cj могут быть нулевыми. В случае, когда C совпадает с C1 или C2, т.е. имеются нижняя или верхняя треугольные матрицы, решения разностных уравнений находятся независимо для каждой компоненты вектора ^ либо сверху вниз (вычисляется вначале &1, а затем &2), либо снизу вверх. При C = C3 система уравнений

&i + тасП2Л&2 = /1,

&2 + та(сП1Л&1 + сП2Л&2) = /2 (3.5)

после исключения компоненты &1 из второго уравнения (3.5) приводится к уравнению относительно &2:

[I + тас^Л - тV^^Л] &2 = /2 - тас^Л/ (3.6)

Его решение может быть получено скалярной прогонкой, после чего явно определяется &1. Реализация схемы (3.1) при C = C4 аналогична рассмотренному выше случаю, но исключение проводится для &1.

Для m = 3 количество представлений матрицы C, для которых возможно решение уравнений (3.1) скалярными прогонками, возрастает. Приведем их вид для случая, когда разностное уравнение (3.4) для отдельной компоненты полученное после исключения

других компонент, содержит лишь первые (с^Л) и вторые (с^ЛстпЛ) разностные производные. Эти матрицы равны:

С =

С4 =

с11 О О 1 с11 с12 с13

с21 с22 О , С2 = О с22 с23

с31 с32 с33 ! О О с33

О с12 О\ 1 О О с13

с21 с22 с23 , С5 = О О с23

О с32 О ! с31 с32 с33

Сз

Сб

с11 с12 с13

с21 О О

с31 О О

О О О

с21 О с31

с31 с32 с33

Таким образом, задача построения безусловно устойчивых разностных схем для решения системы гиперболических уравнений сводится к решению двух подзадач:

1. Для системы дифференциальных уравнений найти все классы схем точной факторизации, т.е. удовлетворяющих условиям (1.6), согласно алгоритму, изложенному в разд. 2.

2. Среди схем в дробных шагах (1.5) или (3.1), полученных при расщеплении матрицы В на В1 и В2, найти все схемы, удовлетворяющие свойству скалярной разрешимости.

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

4. Точная факторизация для одномерных уравнений газовой динамики

Рассмотрим уравнения газовой динамики в консервативных переменных

р

дИ

д£ дх

ри

0, Р = р(р,е)

(4.1)

W

ри2 + р

и2

Е=р(е

где И = | ри

Е / У и(Е + р)

В силу однородности уравнения газовой динамики могут быть представлены в эквивалентных формах

дИ

Ж

д W

д И

В

дИ дИ

дх ' д£ дх '

где матрица Якоби В может быть записана в виде

О

В

дх ,

О

и

V

-(3 - 7(3 - к

и2 , , 3. 2

---ади ад +— ки 7«

к 2

/

Здесь к = 7 — 1; ад = с2/к + 7м2/2. Невырожденным преобразованием матрица В может быть приведена к диагональному виду

В = ЛВЬ,

1

где Б = (и, и — с, и + с)

т.

( 1

Я

V

с

2

1

1

\

2 с2 2с2

и и 1 и 1

С2 2С2 + 2С 2С2 — 2С

и2 и2 1 и и2 1 и

2С2 4С2 + 2к + 2С 4С2 + 2к — 2с/

ь

^ —с2 + к—ки к ^ 2

—си + к

и

2

с — ки к

и2

V

си + к — —с — ки к

/

В соответствии с разд. 2 (см. (2.6) и (2.8)) существуют шесть двухпараметрических расщеплений матрицы Б на Б и Б2 и, следовательно, расщеплений матрицы Б на В?! и В2, удовлетворяющих условиям (1.6). Перебором всех расщеплений Б из разд. 2 можно убедиться в том, что в консервативных переменных р,ри,Е не существует схем с полной факторизацией, реализуемых скалярными прогонками.

Рассмотрим уравнения (4.1) в недивергентной форме

дх

0.

(4.2)

_ ~ дИ ~

Здесь ! — вектор искомых функций; В = А !ВА, где А = -—, В? =

д W

_ Найдем газоди-

д I д И

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

р \ / и р 0

I = | и , В = 0 и 1/р р / \ 0 7р и

и преобразование (4.2) приводит матрицу В к диагональной форме, где

1

и 0 0 Б = ( 0 и — с 0

0 0 и + с

Я

( 0 р

2с 2с2 1

1 0 —-

V

с2

0 -Р —

2с 2с2 )

1

ь

1

сс

— 0 -

рр

с2 0 с2

7Р/р.

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

4.1. Переменные плотность, импульс и давление

Для этих переменных исходные уравнения газовой динамики удобно представить в дивергентно-недивергентной форме

| + В1 = °. (43)

2

с

где

р

ри Р

В

О

д

— и

дх

д

д_

дх

2

дх д

О д_

дх д

V

7Р^

—ис — с —и и— , дх дх дх /

т.е. уравнения неразрывности и движения записаны в дивергентной форме, а уравнение энергии — в недивергентной.

Для уравнений с постоянными коэффициентами матрица В равна

В

О 1 О

-и2 2и 1

22 -ис2 с2 и

а соответствующие матрицы Я и Ь равны

Я

1 1 1

с с

О

р р

с2 О с2 )

Ь

О-

р

1

V

1О О

2с

2с 2с2 1

1

2с2 )

Среди всех расщеплений В согласно разд. 2 существует однопараметрическое семейство расщеплений В на В1 и В2 в виде

/000\ /0 1 0 \

В1

—с/ О

/

В2

-ис2 О и

/с — и2 2и 1 —

/

О

с О

для которых удовлетворяется свойство скалярной разрешимости, а В1В2 = 0 и В = В1 + В2. Здесь / — параметр. Разностная схема

£«,+ 1 _ £«

I + таВ« I + таВ!

—В «£«

(4.4)

для уравнений с постоянными коэффициентами совпадает с нефакторизованной схемой (1.2). Схема в дробных шагах

I + таВ«) Г+1/2

(4.5)

_ _В«£«

= 1«,

; I + таВ«] Г+1 = Г+1/2

£ «+1 = £« + т £«+1

эквивалентна (4.4) после исключения вспомогательного вектора Покажем, что она ре ализуется скалярными прогонками. Пусть

/О 0 0 \ (

вВ« =

ОО

—Л(с/)« О Л с« \ — (ис2)«Л 0 и«Л )

В2«

О

Л

V

—Л(/с — и2)« 2Ли« Л — -(с2 )«Л

О

Г «

О

О

(4.6)

2

2

с

2

2

с

— разностная аппроксимация матричных В и В2 операторов с порядком 0(Л,к), где Л — аппроксимация первых производных, например несимметричным оператором Л = Л- при и > 0 и Л = Л+ при и < 0. Здесь Л± = ±(1 — 7±)/Л,, Т± — оператор сдвига, т.е. Т±/ = Л±ъ а Л — сопряженный с Л оператор. На первом дробном шаге схемы (4.5) разностные уравнения

¿«.+1/2 _ ¿га

¿р ¿р ,

С

та+1/2

m

С + таЛ /

1

сгалга+1/2__¿n+1/2

¿п+1/2 + таигаЛ£п+1/2 = ¿р^ + га(ис2)гаЛ^п+1/2

¿га+1/2 ^га+1/2

разрешаются скалярной прогонкой для ¿р из последнего уравнения и явно для ¿т . На втором дробном шаге система разностных уравнений

с+1 = с+1/2 — ™лс+\

¿m+1 + 2таЛипет+1 + таЛ

(/c - u2)xn+1 + (1 - - ГС+1

С n ¿m,

¿n+1 = ¿n+1/2 - тас2ЛСт+1

решается одной скалярной прогонкой для ¿m+1 из уравнения

[I + 2таЛип - т2а2Л(с2 - и2)пЛ] £

¿n+1 ¿m

¿m+1/2 - таЛ

(/С - u2)ncn+1/2 + (1 --)cn+1/2

(4.7)

полученного исключением ¿n+1 и ¿pn+1 из разностного аналога уравнения движения. После нахождения ¿m^1 явно вычисляются ¿n+1 и ¿p^1.

Таким образом, решение уравнений газовой динамики сводится к двум скалярным прогонкам и явному вычислению остальных функций. Заметим, что для решения уравнений (4.7) при M > 1, т.е. |u| > c, прогонка будет плохо обусловленной и для ее реализации необходимо использовать алгоритм немонотонной прогонки (см., например, [8]). Параметр / может служить для минимизации числа арифметических операций при реализации схемы (4.5). Так, в частности, при / = 0 число арифметических операций на 20% меньше, чем при / = 0. Частный случай расщепления (4.6) при / = 0 рассмотрен в работе [9]. Для уравнений (4.3) с переменными коэффициентами матрица В1В2 содержит один ненулевой коэффициент вида -т2а2и [(с2Л - Лс2)Л)] , т.е. схема (4.5) является схемой с минимальной диссипацией (см. [6]) и схемой точной факторизации лишь при c2 = const.

4.2. Переменные плотность, импульс и поток импульса

Уравнения газовой динамики (4.3) в этих переменных имеют матрицу В следующего вида:

В

0

д_

ôx 0

д д u(u2 - c2)^ (с2 - 3u2)^ dx dx

0 д_

ôx 3 д

dx

Р

pu p + pu2

Р

m R

(4.8)

0

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

О 1 О

В

О

О

1

и(и2 — с2) с2 — 3и2 3и

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

В1

О О

О О

О О

и(и2 — с2) (с — и)(с + 2и) и — с

В2

О О

1 О

0

1

О —и(с + и) 2и + с

О 0 0

В1 = О О О

и(и2 — с2) —2и2 и

В

2=

О О

1 О

0

1

О с2 - и2 2и

В1

О О

О О

О О

и(и2 — с2) (с + и)(с — 2и) и + с

О

1

О

В2 = О О 1

О и(с - и) 2и - с

В1

В1

О О

1 О

1

(с - и) О

В2

\ и(и2 — с2) с2 — 3и2 с + 2и )

(

В1

О О

1\

и

О

О

В2

\ и(и2 — с2) с2 — 3и2 2и /

О О

1 О

1

\

(с + и) О

В2

\ и(и2 — с2) с2 — 3и2 2и — с

(

ОО ОО

1

(

Замечание: В переменных р,ри,р + ри2 разностная схема

£«+1 _ £«

I + таВга )---= —■В«£«

( и - с) 1

\

О О и - с

0 0

и

0 0 1

1 \

\ 0 0 и /

ОО ОО

\ 0 0 и + с )

1

(и + с) 1

(4.9)

может быть реализована скалярной прогонкой без введения расщепления матрицы В. Действительно, система уравнений (4.9), записанных в скалярной форме

+ таЛ£т = —Лт« = —Лт« = /„, + таЛ£д = —Л Я« = /т,

1

£д + 3тампЛ£д + та [(с2 - 3м)пЛ£т + мп(м2 - с2)пЛ^] =

-3мЛЯп - мп(м2 - с2)пЛрп - (с2 - 3и2)пЛт = /д

где = / - /П, сводится (после исключения и £д из уравнения движения) к уравнению

{/ + 3таЛ - т2а2(с2 - 3м2)пЛЛ + т3а3мп(м2 - с2)пЛЛЛ} £д =

:2)пУ 22

= /д - там"(и2 - с2)пЛ / - т2а2Л/т] - та(с2 - 3м2)Л/т,

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

4.3. Переменные плотность, импульс, скорость звука

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

(

В

О

с2

--м

7

2

V

2р

1 О

2м 2ср

(7 - 1)с 7

2р м )

и ее расщепление представляется в виде (/ — параметр)

В1

О

(7 - 1)С

7

(7 - 1)см

V

2р

О

О \

2р1

7

м /

В2

О

с2 - 7м2 + (7 - 1)с/

V

7 О

1

2м

(7 - 1)с 2р

2р

7

(с - 1)

О

/

Матрицы В1 и В2 подобны матрицам для переменных плотность, импульс и давление, т.е. реализация схемы (4.5) идентична рассмотренному выше случаю.

О

4.4. Схема типа предиктор-корректор

Наряду со схемами точной факторизации (1.4) или (1.5) могут быть построены схемы типа предиктор-корректор, где на этапе предиктора исходные уравнения аппроксимируются по схеме расщепления в недивергентной форме, а на этапе корректора восстанавливается консервативность схемы. Такой подход может оказаться более предпочтительным, особенно при расчете разрывных течений. Итак, рассмотрим систему уравнений гиперболического типа, записанную в консервативной (4.1) и неконсервативной (4.2) форме. Разностная схема предиктор-корректор

гга+1/4 _ гп

I-l + В1г+1/4 = о,

та

гга+1/2 — гга+1/4

— + £2Р+1/2 = 0, (4.10)

та

ип+1 — ип

+ ЛЖга+1/2 = 0

т

подобно (4.4), (4.5) аппроксимирует исходные уравнения (4.1) с порядком 0(тт + кк), где т =2 при а = 0.5 + 0(т). Для различного выбора вектора искомых функций ¥ разностные матричные операторы В3 = 1, 2) выбираются из условия скалярной разрешимости схемы на дробных шагах и удовлетворения условий точной факторизации (1.5). Виды таких расщеплений рассмотрены выше в разд. 4.1 — 4.3. После нахождения значений ¥п+1/2, а следовательно, и Wra+1/2 новые значения функций ип+1 определяются явно из последнего уравнения схемы (4.10), аппроксимирующей исходные уравнения в консервативной форме. По построению при а > 0.5 разностная схема (4.10) безусловно устойчива и экономична, так как ее решение сводится к скалярным прогонкам или к схеме бегущего счета подобно рассмотренным выше схемам.

5. Схемы для многомерных задач

Рассмотрим обобщение изложенных выше алгоритмов на многомерный случай. Пусть

д и Л ^ + £

¿=1 3

0

(5.1.)

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

д ¥

Я* N

— + V В

д* + ^ 3 дх

3 = 1 и

0,

(5.2)

полученной из (5.1) при разрешении ее относительно вектора ¥, матрицы В3 имеют действительные собственные числа. Здесь И = И(¥), А = дИ/д¥, W3• = W3•(¥), В3 = А-153А, В, = дW3• /дИ. В частном случае вектор ¥ может совпадать с И. Тогда и каждая одномерная система для ] = 1,..., N

дИ дWj „

или

д* + дх3

— + Вд3

д* дх3

0

(5.3)

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

В

В31 + В

32

(5.4)

— расщепление матрицы В3 на сумму матриц, удовлетворяющую условию точной факторизации (1.6) и свойству скалярной разрешимости разностных схем (см. разд. 3). Тогда факторизованная схема

Г+1 - Г

(/ + таВ}1 Л,-) (/ + таВ^Л,) 7 т 7 + В,Л,-Г? = О (5.5)

или эквивалентная ей схема в дробных шагах типа (4.5) и схема предиктор-корректор вида (4.10) аппроксимируют одномерную систему уравнений (5.3) с порядком 0(тт + Л,), по построению реализуются скалярными прогонками на дробных шагах, безусловно устойчивы при а > 0.5 и удовлетворяют свойству точной факторизации. Заметим, что для каждой одномерной системы уравнений могут быть выбраны свои независимые переменные, т.е. Г, может принимать различные значения при различных направлениях х,. Схемы (4.5) или (5.5) реализуются на временном интервале [п + - 1)/Ж, п + ], где ] = 1,..., N.

Для многомерных уравнений (5.1), (5.2) разностная схема предиктор-корректор с учетом (5.4) может быть представлена в виде

г-1 + В-Л^Г1^ = 0, г-^-+ Вп2Л1^п+1/м = 0

та 11 1 та 12 1

(5.6)

,<п+(М-1)/2М .рп+(М-2)/2М гп+1/2 гп+(М-1)/2М

%_- 1 + Вп Л Гп+(М-1)/2М = 0 1 ' - 1 ' '7 + Вп Л гп+1/2 = 0

--г ВМ1ЛМ% = °--+ ВМ2ЛМ% = 0,

та та

т тп+1 ттп N

+ 5>^'+1/2 = 0.

,=1

Она аппроксимирует исходные уравнения (5.1) с порядком 0(тт + ), по построению на каждом дробном шаге реализуется скалярными прогонками, а на этапе корректора

д

явно. Здесь Л = тах(Л,), Л, = —--+ 0(Лк), т = 2 при а = 0.5 + 0(т), а векторы Г, могут

дх,

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

N гп+1 _ гп

П(1 + таВПА) (/ + таВ^Л,)-т-= -Wn (5.7)

,=1

или эквивалентная ей схема в дробных шагах

= -Wn,

(I + таВ"1Л1) |п+1/2М = |п, (I + таВ"2Л1) |п+1/м = |п+1/2М,

(I + таВп 1ЛМ) -1)/2М = Г+(М-1)/М, (I + таВп2Лм) |п+1/2 = Г+(2М-1)/^,

= г п + т ^п+1

В качестве вектора правых частей выбирается аппроксимация стационарной части урав-

N д

нений (5.2) в недивергентной форме с порядком 0(Л5), т.е. Wn = У2 ВпЛГп, Л? = —--+

с/ с/ с/

,=1 дх,

„п^ „ ^ п д¥ дWj

С учетом соотношения > В7 —— = А > —— аппроксимацию правых частей

3 = 1 дх3 3=1 дх3

N

можно выбирать в предельно дивергентной форме Wn = (А-1)п ^ Л3W3^ Тогда при по-

3=1

лучении стационарного решения методом установления в стационарном случае, т.е. при ¥п+1 = разностная схема (5.7) аппроксимирует стационарные уравнения в консерва-

N N

тивной форме. (А-1) Л3W3l = 0 или ^ Л3W3!' = 0, так как А — невырожденная матри-

3=1 3=1

ца. Обычно значение в выбирается большим (или равным) к с целью повышения порядка аппроксимации. Тогда разностная схема (5.6) аппроксимирует исходные уравнения с порядком 0(тт + тЛк + Л3), а стационарные уравнения после установления — с порядком 0(Л3). Конечно, и в схеме (5.6) на этапе корректора может быть выбрана аппроксимация уравнений в консервативной форме с порядком 0(Л3). Разностные схемы (5.6) или (5.7), как и все факторизованные схемы, безусловно устойчивы при а > 0.5 лишь для двумерного случая и условно устойчивы для N = 3. Заметим, что при в > к устойчивость схемы понижается. Конечно, в многомерном случае схемы (5.6) и (5.7) не будут схемами с точной факторизацией, но они совпадают со схемами типа Бима — Уорминга [10] или Макдональда — Брили [11], которые в отличие от них реализуются векторными прогонками. Изложенная выше технология построения экономичных разностных схем может быть применена при построении схем расщепления на основе метода конечных объемов подобно [12], где предложены схемы приближенной факторизации, обладающие свойством минимальной диссипации, для численного решения уравнений газовой динамики.

Список литературы

[1] Fletcher C.A. J. Computational Techniques for Fluid Dynamics. Vol. I, II. Berlin: Springer-Verl., 1988.

[2] Андерсон Д., Таннехил Дж., Флетчер Р. Вычислительная гидродинамика и теплообмен. Т. 1, 2. М: Мир, 1991. 123 с.

[3] Proc. of 1th — 7th Intern. Symp. on Computational Fluid Dynamics. 1985, 1987, 1991, 1993, 1995, 1997.

[4] Марчук Г.И. Методы расщеплений и переменных направлений. М.: ОВМ АН СССР, 1986. 334 с.

[5] Ковеня В.М., Яненко Н.Н. Метод расщепления в задачах газовой динамики, Новосибирск: Наука, 1981. 304 с.

[6] Ковеня В.М., Лебедев А.С. Модификация метода расщепления для конструирования экономичных разностных схем // Журн. вычисл. математаки и мат. физики. 1994. Т. 34, № 6. С. 886-897.

[7] Ковеня В.М. Методы вычислений (дополнительные главы): Учеб. пособие. Новосибирск: НГУ, 1995. 92 с.

[8] Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978. 501 с.

[9] КовЕня В.М., Грибанов С.В. Об одной схеме расщепления для численного решения уравнений газовой динамики // Вычисл. технологии: Сб. науч. тр. / РАН. Сиб. отд-ние. ИВТ. Т. 2, № 5. 1993. C. 87-94.

[10] BEAM R.M., Warming R.F. An implicit finite-difference algorithm for hyperbolic systems in conservation law form // J. of Comput. Phys. 1976. Vol. 22. P. 87-108.

[11] Briley W.R., McDonald H. Solution of the 3D Compressible Navier — Stokes Equations by an Implicit Technique. Lect Notes in Phys., 1975. Vol. 35.

[12] Ковеня В.М. Схемы расщепления в методе конечных объемов // Журн. вычисл. математики и мат. физики. 2001. Т. 41, № 1. C. 100-113.

Поступила в редакцию 26 февраля 2003 г., в переработанном виде — 6 марта 2003 г.