Оптимальные алгоритмы расщепления для численного решения уравнений Эйлера и Навье - Стокса Текст научной статьи по специальности «Математика»

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

Том 19, № 4, 2014

Оптимальные алгоритмы расщепления для численного решения уравнений Эйлера и Навье — Стокса*

В. М. Ковеня

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

e-mail: kovenya@ict.nsc.ru

Ковеня В.М. Оптимальные алгоритмы расщепления для численного решения уравнений Эйлера и Навье — Стокса // Вычислительные технологии. 2014. Т. 19, № 4. С. 42-60.

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

Ключевые слова: уравнения Эйлера и Навье — Стокса, разностная схема, методы расщепления и факторизации.

Kovenya V.M. Optimum of splitting algorithms for the numerical solutions of Euler and Navier — Stokes equations // Computational Technologies. 2014. Vol. 19, No. 4. P. 42-60.

Numerical solution of Euler and Navier —Stokes equations for compressed heat-conducting gas is considered. We present a class of the implicit differential schemes based on optimum splitting of the initial equations, which is common for the equations written in divergent and not divergent forms. The initial equations are represented in the form of special splitting on physical processes and the spatial directions which influence is minimum. The method is applied for economic differential schemes of predictor type. At the predictor stage, the differential equations are solved on fractional steps, using effective algorithms, for example, scalar sweeps. At the corrector stage, the conservatism of algorithm is restored. The obtained differential schemes approximate the equations with the second order on all variables. They are stable (certainly they are absolutely stable for two-dimensional problems for the appropriate choice of weight parameter and are conditionally steady for spatial problem) thus they are suitable for solution of stationary and non-stationary problems. Their realization is reduced to scalar sweeps unlike implicit differential schemes with directions splitting, which are realized by vector sweeps, where dimension is defined by the equations on space. The analysis of properties of the offered algorithms for the equations of various dimensions is carried out. Offered algorithms allow to choose various gasdynamic variables on fractional steps at the predictor stage that can lead to simplification of their realization and minimize number of arithmetic operations on grid points.

Keywords: Euler and Navier — Stokes' equations, differential scheme, splitting and factorization methods.

*Работа выполнена при поддержке РФФИ (грант № 14-01-00191) и Интеграционных проектов СО РАН № 76 и 130.

Введение

Физико-математические модели, описываемые уравнениями Эйлера и Навье — Стокса сжимаемого теплопроводного газа, используются при численном моделировании в аэродинамике, энергетике, охране окружающей среды, проектировании технических устройств и моделировании технологических процессов. Для их численного решения широко применяются методы конечных разностей и конечных объёмов как наиболее универсальные, пригодные для решения различных классов уравнений (см., например, [1-14]). Увеличение размерности задач и усложнение геометрий расчётных областей накладывают дополнительные требования на численные алгоритмы, которые должны обладать необходимой точностью, иметь достаточный запас устойчивости, удовлетворять свойствам консервативности и быть экономичными с возможностью нахождения решения задачи за разумное время на существующих ЭВМ. Применение явных разностных схем для решения уравнений Навье — Стокса может оказаться неэффективным в силу жёстких ограничений на их устойчивость, особенно при нахождении стационарного решения методом установления, поэтому в данном случае наиболее часто используются неявные разностные схемы, свободные от ограничений на устойчивость или имеющие более слабые ограничения. Для построения неявных разностных схем применяются методы факторизации и расщепления, позволяющие свести решение исходных многомерных задач к последовательному (или параллельному) решению их одномерных аналогов или задач более простой структуры (см. [5, 8-14]). Однако введение расщепления или факторизации приводит к появлению в разностных схемах дополнительных членов второго и более высокого порядка, что может вызвать ухудшение свойств численного алгоритма. Поэтому в задачах численного моделирования проблема построения экономичных алгоритмов остается одной из важнейших.

В настоящей работе экономичные (по числу операций на узел сетки) численные алгоритмы решения уравнений Эйлера и Навье — Стокса строятся на основе метода расщепления, идеология построения которого изложена в [12, 13]. При этом расщепление исходных уравнений выбиралось таким образом, чтобы построенные на его основе алгоритмы были экономичны, обладали большим запасом устойчивости, а влияние расщепления было минимальным.

1. Алгоритмы расщепления для уравнений газовой динамики

Рассмотрим систему уравнений Эйлера при отсутствии внешних сил (см., например,

[12, 14])

ди Л дЖ3 . ,

+ ^ 1ПГз = °, Р = Р(Р,е),

дЬ

дхз з=1 з

где хз — декартовы координаты,

и

( Р \

рУ] РУ2 рУз

V Е )

W,

( РУ3 \ РУ1 Уз + 5] р ру2Уз + 5]р

рузуз + 63Р

V Уз (Е + Р)

, Е = р(е + у7 2),

при размерности задачи по простраству N = 3. Далее для определённости уравнение состояния будем задавать в виде р = (7 — 1)ре. Пусть Л7 = д/дж7- + 0(Л,к) — аппроксимация производных с порядком к. Операторы Л7 могут задаваться на симметричных или несимметричных шаблонах [3, 12]. Разностная схема с весами

ип+1 _ туп N

-+ Л,- [а¥"+1 + вWn] =0, в =1 — а,

.7 = 1

аппроксимирует уравнения (1) с порядком 0(тт + Нк), где при а = 0.5 + 0(т) т =2, иначе — т =1. Эта схема нелинейна при а = 0 и требует итераций. Линеаризуем векторы и, W относительно вектора С по формулам

Ип+1 = Ип + т (дИ/дС)п (дf /д£)п + 0(т2) = Ип + Лп(:Т+1 — Г) + ..., А = дИ/дС,

Wn+1 = wn + т (дWj/дС)п (дС/д£)п + 0(т2) = Wn + Бп(:Т+1 — Г) +..., В, = дWj/д£. Тогда разностная схема

сп+1 _ сп М N

С-+ ^ Л, Wn = 0, С = Лп + та£ Л, Бп

т 7=1 7=1

является линейной, аппроксимирует уравнения (1) с порядком 0(т + кк) и реализуется матричными прогонками. Для получения экономичных схем обычно применяются методы расщепления по направлениям или факторизации (см. [3-5, 9-14]), полученные

при замене стабилизирующего оператора С на сумму одномерных операторов

М 1

С = Е( ^ ЛП + таЛ7 7

7=1

Разностная схема расщепления

сп+1/М _ с п

(Лп + таЛда-+ Л^+1/М = 0,

т

сп+1 _ сп+(М-1)/М

(Ап + таЛм БN)-+ Лм WN+1 = 0

также аппроксимирует уравнения (1) с порядком 0(т + Л,к), но на каждом дробном шаге реализуется уже векторными прогонками. При а > 0.5 она безусловно устойчива для уравнений любой размерности. Напомним, что схемы расщепления при N > 2 не являются схемами полной аппроксимации, так как не приводятся к канонической форме [8, 14].

Схема приближённой факторизации [9, 10]

сп+1 _ сп N

(Лп + таЛ1Вп) ■ (А-1)п■ (Лп + таЛ2Вп) ■ ■ ■ (Лп + таЛмBN)-+ ^ Л.W7n = 0, (3)

т 7=1

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

N

АП£П = _ ^ л,,

3 = 1

(Ап + таЛ^)^1^ = Ап|п, (Ап + тaЛN BN )Г+1 = Ап|п+^-1)/N,

£П+1 _ £П + ^П+1

аппроксимирует уравнения газовой динамики с тем же порядком, что и схема (2), и также реализуется векторными прогонками, но она является схемой с полной аппроксимацией [3, 8], безусловно устойчива при а > 0.5 для двумерных уравнений и, как всякая схема приближённой факторизации для гиперболических уравнений, условно устойчива в трёхмерном случае [12]. Разностные схемы могут быть построены и для решения уравнений газовой динамики, записанных в недивергентной форме

N

I ^ в, г

0, в,

,=1

А-1 д^з А М

д дхз,

А

с и

Ж,

разрешённых относительно вектора ¥, где его выбор определяет форму матричных операторов В,. Наиболее простую форму уравнения газовой динамики (1) принимают в переменных плотность р, скорость V и давление р, т.е. для вектора ¥ = (р,у1,у2,у3,р)т. Тогда операторы В, равны

(с

в,

( у, б)р б,3р 0 \

0 у, 0 Ь,1 д

0 0 у, Ь.,3 дх,

V 0 с1 с, с3 с, у

или

в,

схзУз

V

0 0

0

_д_

дх,

1 д

С1 -

3 дхз

0 0

д

3 дх,

А

, дх,

б1 —

, дх

0

д

Ь3

,

д_

1 дх, д,

у

дх, /

если уравнение неразрывности записать в дивергентной форме. Здесь Ь, = б,/р, с, = б,7р. Аппроксимируем операторы В, сеточными операторами В,^ с порядком 0(кк) по формулам

( Л, у, 0

0 0

с1Л

0 0

0.

Ь1Л

с3Л, у,Л, )

Здесь, как и в [12], первые производные аппроксимированы операторами Л, = д/дх ■ + 0(^) с учётом знака скорости у, с первым порядком аппроксимации

Л, = Л,т, Л, = Л,±, если у, > 0, и Л, = Л,

Л, = если у, < 0,

0

V

0

0

0

или симметричными операторами Л7- = Л7 = (Л1- + Л1+)/2, где Л7^ = ±(1 — Т±)/Л7-а = /±1 — операторы сдвига в направлении ж7-.

Разностная схема приближенной факторизации

N сп+1 _ сп N

П [! + та в;л]-= — ^ Б^Г (5

7=1 т 7=1

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

N

Г = — У! Б"ЛСп,

7=1

(I + та Б^)^1^ = Г ... (I + та Б^ )|п+1 = -1)/N, (6)

£п+1 = £ п + т ^п+1

где ^ = (СР,С1,С2,Сз,Ср)т — вектор невязок, аппроксимирует уравнения (4) с порядком 0 (т + Лк), так как коэффициенты матриц Бп аппроксимированы на нижнем временном слое п, и также реализуется на дробных шагах векторными прогонками. При решении стационарных задач методом установления более предпочтительной может оказаться разностная схема

N сп+1 _ сп N

Д[1 + та Б*Л]-= — (Лп)-1 ^ Л7 w7n,

7=1 т 7=1

полученная из (5) при замене оператора правой части в недивергентной форме на оператор в предельно-дивергентной форме, как это следует из соотношений

® _ —Л-1 V _ — у Б с

дг = Л 917 = 2-,Б7 ' ■

7=1 7=1

Для повышения порядка аппроксимации схемы (5) разностные операторы Бк^ на дробных шагах и в операторе правой части могут задаваться с различным порядком. Такую схему назовем схемой с несогласованной аппроксимацией. В отличие от (5) она аппроксимирует исходные уравнения с порядком 0(т + тЛк + Л1) и реализуется при к =1, I = 2 на дробных шагах трёхточечными векторными прогонками. При установлении она аппроксимирует стационарные уравнения

N

£ Л7 = 0

7=1

в консервативной форме с порядком 0(Л1). Аналогично строится схема расщепления

сп+1^ _ сп

(I + таЛ^)с-- + (Л-1)nЛlWn+1/N = 0,

1 т 1

сп+1 _ с-1)^

(I + таЛN БN)-+ (Л-1)^ WN+1 = 0.

т

Заметим, что реализация разностных схем при переходе с п-го слоя на п + 1-й слой требует Мш3 операций на узел сетки, где ш — количество уравнений.

Для получения экономичных схем, реализуемых скалярными прогонками, в работе [12] было введено расщепление уравнений по пространственным направлениям и физическим процессам. Исходя из характера гидродинамических сил, векторы потоков подразделялись на потоки, обусловленные градиентами давления и действием конвективных сил по каждому пространственному направлению, т.е. исходные уравнения (1) представлялись в виде

д и ^ ^ д

дг

з=1 1=1

дхз

где

( РЩз \

ру2Уз рМзЪз

\ Е )

( 0 \

8)р

V ЩзР )

что позволило исходную многомерную задачу представить в виде совокупности более простых одномерных задач, описывающих элементарные физические процессы:

1 д И д—3

2М дг

дх3

0 = 1,. . . ,М, I = 1, 2, N = 3),

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

| + ££в; Г = 0, (8)

3=1 1=1

где операторы В3 для £ = (р, щ1, щ2, щ3,р)Т равны

В]

00 00

00

V0 с]

00

0 Ь]

0 Ь3

3 0 )

д дхз з

В2

( д

дхз з 0

0 0

'дхз

0 0

0 0

дхз

0 0

д дхз

Как и выше, аппроксимируем операторы Вз сеточными операторами В^ по формула:

м

0 0 0 0 .0 .0 0 Ь]Лз

0 0 0 сз Лз .0 с3Л • сз Лз Ь3Лз 0

В%

( Лз V 0

0 0

0 .0 0 \

Щ Лз .0 0

0 Щз Лз 0

0 .0 Щз Лз )

0

0

V

0

V

0

V

Тогда разностная схема приближенной факторизации

N сп+1 _ сп N

Па + та Б^а + та Б^)-т-= — (Л-1)п £ Л7(10

7=1 т 7=1

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

= —(Л-1)п £ Л7 W п,

N

п £ ^wп,

7=1

(I + таБ1Лп+1/^ = Г, (I + таБ^)^1^ = ,

(I + таБ^)Г+(^-1)/2N = -1)/N, (I + таБ^)Г+1 = -1)/2N,

сп+1 = сп + т|п+1 (11)

аппроксимирует исходные уравнения (7) или (8) с порядком 0(т+ Лк), но реализуется на дробных шагах скалярными прогонками. Действительно, на нечётных шагах решение системы разностных уравнений (для ] = 1,..., N)

.п+(27-1)/2^ _ п+(7-1)^ СР = ?Р ,

еп+(27-1)/^ + таЬ1Л7срп+(27-1)/^ = еп+(7-1)^ (I = 1, 2, 3),

£Рп+(27-1)/27 + тас7 Л7 еп+(27-1)/™ = ерп+(7-1)/2

„ n+(27-1)/2N

после исключения компонент скоростей V из последнего уравнения сводится

к уравнению

[I — т2а2с7Л7а7Л7] Срп+(27-1)/^ = Срп+(7-1)/2 — тас7Л7Сп+(7-1)^,

решение которого находится скалярными прогонками (по каждому пространственному направлению), после чего явно вычисляются невязки скоростей, а значение плотности переносится с предыдущего шага.

На чётных шагах разностные уравнения

[1 + таЛ7 ^ Ср7 = СР

„п+7/^ _ ¿п+(27-1)^ 5Р = СР

[1 + та^] С = СГ+(27-1)/2^,

[1 + та^-] СГ7^ = СРп+(27-1)/^ (I = 1, 2, 3)

ъР

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

сп+1/2^ _ сп сп+1/^ _ fn+1/2N

с-- + Б^1^ = 0, с-1-+ Б?. = 0,

та та

сп+^-1)/2^ _ fn+(N-1)/N

с-Г-+ -1)/^ = 0,

та

сп+1/2 _ fn+(2N-1)/2N

с-с-+ Б^^2 = 0. (12)

та

Эта схема аппроксимирует исходные уравнения с порядком 0(т + Нк) и реализуется подобно (11) скалярными прогонками. Заметим, что число прогонок в схемах (10), (12) с расщеплением операторов по физическим процессам и пространственным переменным (9) равно N(т + 1), а число операций на узел расчётной сетки при переходе с п-го на п + 1-й слой составляет Я = N (т +1) в отличие от схем с расщеплением операторов по пространственным направлениям, где оно равно Я = Nm3.

Заметим, что введённое расщепление (9) для уравнений в недивергентном виде (8) не эквивалентно их расщеплению в дивергентном виде (7). Действительно, уравнениям (7) соответствует недивергентная форма расщепления в виде

в]

( 0

о

о

V0

о 0

о

С2 )

0 0

о

сз С2 )

о \

Ь1

3

Ь3 о

(

д дх);

в2

vз 5)р 5)р о

о v) о о

о о vз о

о с1) с3 с1) v)

дх3

где

в!

дWlj

А-1_3,

дх3

Ь) = 5)/р, с1з = 5)р, ¿2з = (7 - 1)3

При данном расщеплении операторы В) в отличие от (9) содержат как конвективные члены уравнений, так и часть членов с давлением. Поэтому разностные схемы (10) или (12) на дробных чётных шагах реализуются векторными прогонками. Подобная ситуация (неэквивалентность расщепления в дивергентном и недивергентном виде) возникает и при выборе других искомых функций. Как следствие, это приводит к различным разностным схемам, свойства которых также различны, что проявилось уже при решении одномерных задач (см. [16]). Остановимся на данном обстоятельстве более подробно и рассмотрим, каким расщеплениям уравнений Эйлера в дивергентной форме должны соответствовать их расщепления в недивергентном виде. При этом рассмотрим те из них, при которых разностные схемы являются безусловно устойчивыми (или имеют слабые ограничения на устойчивость), реализуются на дробных шагах скалярными прогонками, а влияние расщепления на свойства разностных схем по сравнению со схемами расщепления по направлениям минимально.

2. Оптимальные алгоритмы для уравнений газовой динамики

Представим систему уравнений газовой динамики (1) в виде расщепления по физическим процессам и пространственным переменным в дивергентной

д Wj

ЯТТ N 2

д и ^ ^

~Ж + ^^ ~дх 3=1 1=1

о, где W]

или недивергентной форме

N

дЬ

+ ££ в) £

)=] 1=1

(

о

531р

5332р

533р

\

w2

\ V) (ре + р) /

о, где В) = А

1

Pvз

рvlv) р^)

рш^ V рvзV2/2 )

д W) д д£ дх)'

:13)

Вектор учитывает силы, обусловленные градиентом давления, а '2 — конвективные силы (сравнить с расщеплением в (7)). Покажем, что именно данное расщепление соответствует принятым выше требованиям к разностным схемам. Заметим, что расщеплённые системы уравнений

1 Ш Д д'2 ---+ ^--

2 дЬ -=1 дх-

1 ди Д д'2

2 дЬ -=1 дх-

0,

где и = (р, ру^^, ру2, ру3,,р)1 , = (ру-,ру1у-, ру2V-,ру3у-, 0)т, эквивалентны, и это позволяет при построении разностных схем рассматривать различные газодинамические

1 дЕ . 1 N д о 0

переменные. Действительно, из уравнения для полной энергии ———+— ^

2 дЬ 2 -=1 дх-

ру- у

следует, что д(ре)/дЬ = 0 или др/дЬ = 0. Поэтому при решении уравнений (13) на различных шагах могут выбираться различные переменные £ = (р,у1,у2,у3,р)т или £ = (р, ру1, ру2, ру3,р)т. Пусть £ = (р,у1,у2,у3,р)т. Тогда операторы В- в (13) равны

В1

00 00

00

V0 с)

0 0

0

ь1

-

0 -

-у-

д дх-'

В2

(

дх- -

V

0 у-

0 0

дх-

0 0

00 00

^ дх- 0 00

\

/

'14)

Расщепление операторов (14) отличается от введённого ранее расщепления (7) лишь для уравнения энергии: конвективный член в этом уравнении из оператора В2 перенесён в В1, т.е. при расщеплении (13) все члены уравнения энергии содержатся только в операторе В]. Как следствие, это приводит к упрощению решения. Пусть В-^ = В- + 0(Нк). Тогда

В^

0 0 0 0 \ / Л-у- 0 0 0

0 0 0 61Л- 0 у-Л- 0 0

, В2 =

0 0 0 Ь3Л- 0 0 у-Л- 0

0 с1Л • с3Л • у-Л- / V 0 0 0 0

Разностные схемы (10) или (12) с расщеплением (14) аппроксимируют систему уравнений (13) с порядком 0(т+кк) и реализуются на дробных шагах скалярными прогонками подобно схемам (10), (12). При этом число прогонок равно Мш (по одной прогонке на нечётных шагах и по ш — 1 прогонкам на чётных шагах по каждому направлению) в отличие от введённого расщепления (9), где оно равняется N(ш + 1). Расщепление уравнений назовем оптимальным, если соответствующие неявные разностные схемы реализуются экономичными скалярными прогонками, причём число таких прогонок (и, как следствие, число арифметических операций) при переходе с одного временного слоя на другой не больше числа уравнений. Этому требованию удовлетворяют схемы (10), (12) с расщеплением (14).

0

0

з

с

з

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

1п+1/42 _ 1п 1п+1/22 _ 1п+1/42

1-- + вЦ+1/42 = 0, 1-1-+ ВЦ+1/22 = 0,

та та

1п+(22-1)/42 _ 1п+(2-1)/22

1-1-+ в21п+(22-1)/42 = 0,

та

^п+1/2 _ £п+(22-1)/42

' 1

та

+ В2ЛЛп+1/2 = 0,

т тп+1 ттп 2

и т- и + £Л= 0, (15)

Т 3=1

где = В, + О (кк), Л,Wjh = дWj/дх, + О (к1). Схема (15) аппроксимирует уравнения (1) с порядком О(тт + ткк + к1), где т = 2, при а = 0 . 5 + О(т), консервативна и, в частности, при аппроксимации операторов на дробных шагах с первым порядком для к = 1, а на этапе корректора со вторым порядком (при I = 2) аппроксимирует исходные уравнения со вторым порядком по всем переменным. Для линеаризованных уравнений газовой динамики схема (15) совпадает со схемой (10), т.е. она безусловно устойчива при а > 0 . 5 в двумерном случае и условно устойчива в трёхмерном. Выбор различных газодинамических переменных на этапе предиктора может приводить к разностной схеме, реализация которой на дробных шагах наиболее проста. Действительно, для вектора 1 = (р,у1,у2,у3,р)Т на нечётных дробных шагах решается система разностных уравнений (для ] = 1, 2, 3)

рп+(23—1)/42 = рп+(3—1)/22,

уп+(2з—1)/42 + таЬЗЛ,рп+(23—1)/42 = ^п+(з—1)/22 (I = 1, 2, 3),

рп+(2з—1)/42 + таЛзрп+(2з—1)/42 + та £ 4Л,-у^-1^42 = рп+(2,-1)/22 .

1=1

-Л- п+(23-1)/42

После исключения Уг из уравнения для давления их решение сводится к ре-

шению уравнения

[I + тау,Л, - т2а2с,Л,ЦЛ,] рп+(2,-1)/42 = рп+(,-1)/22 - тас3Л,уп+(з-1)/22

скалярными прогонками (по каждому направлению). Затем явно вычисляются новые значения скоростей, а значение плотности переносится с предыдущего шага. На чётных шагах

[1 + так,уп\ рп+,/22 = рп+(2,-1)/42, [1 + таупЛ,] уп+,/22 = уп+(2,-1)/42 (I = 1, 2, 3),

рп+,/22 = рп+(2,-1)/42

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

ных

1 = (р, ру1, ру2, ру3,р)Т на чётных шагах разностные уравнения аппроксимируются

в дивергентной форме

[1 + та Л, V?] рп+з/2М = рп+(2,-1)/^ , [1 + та Л,V?] (рт)^'^ = (р^)п+(27-1)/4^ (/ = 1, 2, 3),

рП+]/2И = рП+(2,-1)/4М

и решаются четырьмя скалярными прогонками. На этапе корректора схемы (15) новые

N

значения функций находятся по явным формулам Ип+1 = ип — Л,"^П+1/2. Таким

,=1

образом, решение уравнений газовой динамики при переходе с п-го шага на п + 1-й шаг может быть получено за Ыт скалярных прогонок. Для линеаризованных уравнений разностные схемы (12) и (15) безусловно устойчивы при а > 0.5 для двумерных уравнений и условно устойчивы в трёхмерном случае (см. ниже раздел 4). Аналогично строятся экономичные разностные схемы и для других переменных.

Замечание. Аппроксимация производных в векторах дWз• /дх, операторами Л, Wjh со вторым порядком приводит к немонотонным схемам. Для устранения осцилляций подобно [17, 20] могут быть введены сглаживающие операторы.

3. Оптимальные алгоритмы расщепления для уравнений Навье — Стокса

Рассмотрим систему уравнений Навье — Стокса сжимаемого теплопроводного газа в декартовых координатах х, (см., например, [12, 15]) при отсутствии внешних сил:

N

ди+^ д^ = о,

дг

дх,

:1б)

где

и

р

рщ

рУ2 рУз

V Е I

3 ди ■

Е = р(е + и2/2),а1у V = £

W,

рцз

ри1 V, + 51р — а,1 рУ2 V, + , — а"2 риз V, + б^р — а,

V

V, (Е + р) — к

дТ

дх,

а,

/

з

Т V ,=1

2, а,

з дv ■ дv

!>*а5, а, = л v+у( дх^ + дх,

)

,=1 дх, ,= 1 г=1 дх,

девиатор тензора вязких напряжений, Л = у' — (2/3)у; у, у' — коэффициенты первой и второй вязкости; к — коэффициент теплопроводности. Для замыкания системы уравнений (1б) задаются уравнения состояния

Р = P(P,e), е = е(р,Т),

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

2

Проведём подобно [12] анализ уравнений Навье — Стокса исходя из характера гидродинамических сил. Последние могут быть подразделены на инерционные, или конвективные, силы и силы, обусловленные градиентом давления и диссипативными эффектами вязкости и теплопроводности. В соответствии с этим векторы потоков можно представить в виде суммы векторов

= W31 + W32 +

где

Wi

з

( РУз \

рУгУз

ру2Уз РУ3У3

V УзЕ

Уз и,

W2

з

( 0 \

w3 = -

з

а

\ к(дТ/дхз) + аз

\ Узр )

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

дИ

~дЬ

+

N 2

£ £

з=1 г=1

д W5

дхз

:17)

где по аналогии с уравнениями газовой динамики введено расщепление уравнений по физическим процессам в специальной форме:

(

W1

0

зр

з

\

(

w2

V Уз (ре + р) - к(дТ/дхз) /

рУз

рУгУз - аз1 2

рУ2Уз - а2

\

рУзУз - а'

\ РУзУ72 - аз )

:18)

Здесь аз = £ Утат, ат = [^ + ^ + л)]

дУт дхз

операторы Wj, W2 являются одно-

т=1 дхз

мерными и содержат все газодинамические члены (как для уравнений Эйлера) и часть

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

1 дИ д Wг1

2Ы дЬ

дхз

0 и уравнения

1 дИ

т +1 ~т

= Е. Тогда по аналогии с уравнениями

Эйлера уравнения Навье — Стокса (17) можно записать в форме расщепления по пространственным направлениям и физическим процессам в недивергентном виде

Р>Г N 2

I+£ £ вз£

з=1 1=1

Е, где Вз

А-1 дW3• д

д£ дхз'

Е = А-1Е

Л ди А = д* .

:19)

0

1

Как и выше, в качестве искомых функций будем выбирать плотность, скорость и давление или плотность, скорость и температуру. Тогда для р = (7 — 1)ре и вектора ¥ = (р^,v2,v3,p)T операторы равны

В1

в.2

0 0

0 0

д

0 0

0 0

1 д

С1 -

, дх,

дх,Vз 0

0 0

дх,

— $

0 0

0

»Л

, дх

\

=зА

, дх,

V

,

ъзА

, дх, , дх, ,

0 0

0 0

д

ЧЩ — 0 00

где

5\ тр, й,

1д

д

р дх, , дх,

,1м , л ,4 (7 — 1) д

у + 5,(Л + у), $4

д1

к

си дх, дх, р

Подобный вид операторы в, принимают и в других переменных. Аппроксимируем вторые производные й, разностным операторами со вторым порядком, а операторы В, операторами В,ъ:

0 0

0 0

00

\0 с) Л,

0 0

С3Л-

\

Л, ъ}л

Л,

VЛ, — , )

( Л,V

В22>г

V

,, 0

0 0

0

vзЛ, — ,

0 0

0 0

0

0 0

0 0

(20)

Разностные схемы приближённой факторизации (11) или предиктор-корректор (15) с расщеплением операторов в виде (20) аппроксимируют уравнения Навье — Стокса (16) с порядком 0(т + Ьк) и 0(тт + тЬк + Ь), где т = 2 при а = 0.5 + 0(т). Как следует из вида матричных операторов их структура подобна форме операторов для уравнений газовой динамики. Поэтому разностные схемы (11), (15) для уравнений Навье — Стокса реализуются соответственно Ыт скалярными прогонками. Для уравнений Навье — Стокса более удобным в качестве переменных представляется использовать плотность, скорость и температуру, для которых обычно и задаются краевые условия. Как пример рассмотрим реализацию схемы предиктор-корректор (15) для вектора

0

0

0

0

V

ъ

,

0

0

¥ = (р, v1,v2 , v3, Т)т и уравнения состояния р = (7 — 1)рТ. На нечётных дробных шагах решение системы разностных уравнений

рп+(2,-1)^ = рп+0-1)^

V,

п+(2,-1)/^ + та,Л,(аТ)п+(2,-1)/^ = vn+(з-1)/2N (/ = 1, 2, 3),

1 I

Тп+(2,-1)^ + таЛ,Тп+(2,-1)^ + та с1 Л,Vn+(2з-1)/4N = Тп+(2,-1)^

1=1

где Ъ, = 5,/р, а = (7 — 1)р, сводится, как и при решении уравнений для газовой динамики, к скалярной прогонке при нахождении температуры

[I + таvзЛ, — т2а2с,Л,ЦЛ,а] Тп+(2,-1)/^ = Тп+(,-1)/^ — та^Л,vn+(з-1)/2N,

n+(2j-1)/4N

после чего явно вычисляются значения ^ , а плотность переносится с преды-

дущего шага.

На чётных шагах разностные уравнения

[1 + таЛ,V"] р^^ = рп+(2,-1)/^, [1 + та,, — vГ+з/2N = vГ+(2з-1)/4N (1 = 1, 2, 3),

Т n+j/2N = Т п+(2,- 1)/4N

подобно схемам для уравнений газовой динамики решаются скалярными прогонками независимо для каждого уравнения, кроме уравнения для температуры, значение которой переносится с предыдущего шага. Новые значения функций на слое п + 1 вычисляются явно на этапе корректора схемы (15). Таким образом, рассмотренные разностные схемы являются оптимальными по числу операций при переходе с одного временного шага на другой. Аналогично могут быть построены разностные схемы и для других газодинамических переменных. При несогласованной аппроксимации операторов на этапах предиктора и корректора необходимо, как и для схем для уравнений газовой динамики, вводить сглаживание операторов на этапе корректора. Заметим, что введение сглаживания при решении многомерных уравнений приводит к практическому отсутствию осцилляций, но, как показывают расчёты [17-20], устойчивость схем понижается, хотя расчёты можно проводить и при числах Куранта К > 1.

Замечание. Метод расщепления может быть применён и при построении экономичных алгоритмов для уравнений более общей структуры, описывающих более сложные процессы, например, влияние химических реакций. В этом случае возможны два подхода. Первый подход состоит в рассмотрении расширенной системы уравнений, содержащих наряду с законами сохранения балансов массы, импульсов и полной энергии и законы сохранения для отдельных компонентов смеси. Для расширенной системы уравнений может также вводиться расщепление уравнений по физическим процессам и пространственным направлениям. При этом размерность матричных операторов возрастает с увеличением размерности вектора искомых функций И или ¥. Далее подобно схемам для уравнений Навье — Стокса строятся разностные схемы типа приближённой факторизации или предиктор-корректор. Их устойчивость будет зависеть от вида расщепления. Примеры построения таких схем приведены в [20] для численного решения уравнений магнитной гидродинамики. Второй подход состоит в рассмотрении расширенной системы как совокупности двух подсистем, взаимосвязанных через правые

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

Реализация краевых условий в алгоритмах расщепления представляет собой отдельную проблему. На сегодня нет строгого обоснования реализации краевых условий в неявных схемах расщепления. Их изложение для численного решения уравнений Эйлера или Навье — Стокса не представляется возможным без привязки к конкретной краевой задаче. Некоторые случаи неявной реализации краевых условий для уравнений Эйлера и Навье — Стокса приведены в работах [12, 14, 21].

4. Анализ свойств алгоритмов расщепления

Выше отмечалось, что расщепление по пространственным направлениям приводит к появлению в разностных схемах дополнительных членов, отсутствующих в исходной схеме с весами. Остановимся на вопросе влияния расщепления на свойства разностных схем на примере уравнений газовой динамики. Исследование свойств разностных схем проведём для уравнений (1) в переменных - = (р,у1}у2,У3,р)Т. Разностная схема расщепления по пространственным направлениям (5) после исключения дробных шагов может быть представлена в виде

-га+1 _ £п 2

-+ У В^(аГ+1 + в П = Ф (21)

т *—'

.7 = 1

или

2 -п+1 _ -га 2

П + та В;Л]-т-= - £ В— + Фь (22)

7=1 Т 7=1

где

а2—га+1 _ о2-га

Ф = т2(В^ь + ... + В2-1ьВ*ь)-— + Т2ВШ ■ ... ■ В*ь(а3Г+1 + в-п) + ...,

т

Ф1 = (2а - 1)т(В^В^ + ... + Вм-шВм— + (3а2 - 3а + 1)т2ВШ ■ ... ■ В*— + ...

Погрешность аппроксимации Ф = 0(т2) и есть те добавочные члены, которые возникли за счёт расщепления исходных операторов. Напомним, что все схемы расщепления для трёхмерных задач не приводятся к каноническому виду, так как Ф1 = 0 для всех значений весового параметра а при N > 2. Это является недостатком схем расщепления при решении пространственных задач, особенно при нахождении стационарного решения методом установления, поскольку стационарное решение будет зависеть от шага т. В отличие от них схемы приближённой факторизации и предиктор-корректор являются схемами с полной аппроксимацией. Поэтому их использование является предпочтительным, так как заранее нельзя предположить, будет ли решение стационарным даже при стационарных краевых условиях.

Остановимся на анализе свойств схем приближённой факторизации и проведём сравнение схемы (5) с расщеплением уравнений по пространственным направлениям

N

Д[1 + та в;л]

7=1

—га+1 _ — п

N

- £ в;—«

7=1

(23)

и схемы (10) с расщеплением уравнений по физическим процессам и пространственным направлениям в виде (9) и (14)

С

—п+1 _ — п

N N

£ В^Г, С = П(1 + тав;л)(1 + таВ'Л) .

7=1 7=1

(24)

Напомним, что схема (5) реализуется на дробных шагах векторными прогонками, а (10) -скалярными. Введение расщепления в (24) приводит к появлению дополнительных чле-

2 2т->1 т->2

нов вида т2а2В7^В7^, что следует из соотношений

N

N

N

с = П [! + таВ7-Л + т2а2в1^] = Д [I + таВ^]+фо, Фо = т2а2 £ В^+0(т3).

7=1 7=1 7=1

Оценим их влияние при различной форме расщеплений. Эти дополнительные члены для расщепления (9) и (14) соответственно равны

Фо

Фо

( 0 0 .0 0 \

N 0 51 .0 0

т 2а2 £ в;в; = 7=1 0 0 ^ . со 0

\ 0 ¿3 0 /

0 0 .0 0 \

N 0 0 .0 0

т2 а2 £ В^ =

7=1 0 0 .0 0

V 0 ¿1 ¿3 0 /

где 4 = Ь7Л7-VЛ7, ¿7 = 7Л7VЛ7. Таким образом, при расщеплении (9) дополнительные члены присутствуют во всех уравнениях, кроме уравнения неразрывности. Во втором случае введение расщепления (14) приводит к их появлению лишь в уравнении энергии. Кроме того, разностные схемы с расщеплением операторов в виде (14) при переходе с п-го на п + 1-й шаг реализуются за Ыт скалярных прогонок в отличие от (9), где число прогонок равно N (т +1).

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

N

Г = 7ь..7 = №ей, Л = е-, Я = г^к3 к3,

7=1

где i — мнимое число. Очевидно, fn+1 = Л • fn и ||fn+1|| = |Л| • || max I fn I. Тогда тЛ,-fn = idjfn, где А,- — симметричный оператор, d

Характеристическое уравнение схемы (24) примет вид

, где || fn|| = (т/hj )sin(fcj hj).

det

3

П1 2 X—^ l 2

(I + iadj B^XI + iadj Б^)(Л - 1) + i ¿j dj (B^ + B^) j=i j=i

0.

(25)

Получить решение уравнений (25) в общем случае затруднительно, поэтому рассмотрим предельные случаи решений: у ^ в (можно пренебречь конвективным переносом) и у ^ в, где в = \fjpJP — скорость звука. В первом случае характеристическое уравнение схемы (25) примет вид

Л1 = 1, det

Л - 1

0 Л - 1

00 id^2(a^ + в) id2s2(aЛ + в)

-a2s2d1d2^ - 1) -a2s2d1d3^ - 1) id1(aЛ + в) -a2s2d2d3^ - 1) id2(aЛ + в) Л - 1 + в)

id3s2 (аЛ + в)

Л1

Его корни равны

Л

1,2,3

1, Л

4,5

1 + а(1 - a)R2 ± i(2a - 1)R

1 + a2R2 ,

где Я2 = в2 ^2 + d22 + d2i + а2в2(d2d2 + d21d23 + + а4в4При а > 0.5 получим, что |Л|г < 1(1 = 1,...,5), т.е. разностная схема (24) безусловно устойчива. Во втором предельном случае характеристическое уравнение (25) можно представить в виде

det

П(1 + iadjVj)(Л - 1) + j=1 j=1

j vj

0, Л5 = 1,

где I — единичная матрица 4 х 4. Корни (25) равны

A + i(B - V)

Л

1,...,4

A + iB

Л5

где А =1 — а2(у2+ + В = а(У — а2у2ь2V = ^ djу.

У=2

Пусть для некоторых узлов сетки выполняются условия dj Уу = do = (тУу/Ну) 8т(к Ну). Тогда удовлетворение условия устойчивости схемы |Л2,...,5|2 < 1 приводит к неравенству 3d2 < 2d0B = 2аd2 [3 — (аd0)2]. Это означает, что временной шаг т должен удовлетворять соотношению

т < h3(2a - 1)/2а3у/v0, v0 = max |vj|

и максимальное значение шаг т принимает при значении весового множителя а = 3/4. В данном случае условие устойчивости принимает вид т < 4Н/(3ь0), т. е. близко к устойчивости явных схем. Таким образом, разностная схема приближённой факторизации (24) при преобладании конвективного переноса условно устойчива для всех значений весового параметра а, а при а = 0.5 абсолютно неустойчива. Следовательно, схемы

n

0

I

1

приближённой факторизации и предиктор-корректор для трёхмерных уравнений газовой динамики условно устойчивы. Подобные оценки устойчивости были получены в [12] для схемы типа (24) при расщеплении в форме (9). Отметим, что для двумерных уравнений разностные схемы приближённой факторизации и схема предиктор-корректор безусловно устойчивы при a > 0 . 5.

По предложенным алгоритмам были проведены некоторые расчёты сверхзвуковых вязких течений в канале и около тел с возникновением зон ударного перехода, их взаимодействием между собой и появлением зон отрыва (см. [18-20]), подтвердившие эффективность алгоритмов и достаточную точность расчётов даже без привлечения алгоритмов сгущения сеток. Простота реализации представленных схем, возможность изменения параметров потока и шагов сетки в широком диапазоне на основании большого запаса устойчивости схемы, экономичность алгоритма — всё это позволяет использовать предложенные алгоритмы для численного решения различных классов задач.

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

[1] Роуч П. Вычислительная гидродинамика. М.: Мир, 1980. 616 с.

[2] Флетчер К. Вычислительные методы в динамике жидкостей: В 2-х томах. М.: Мир, 1991. 504 с.

[3] Марчук Г.И. Методы вычислительной математики. М.: Наука, 1980. 536 с.

[4] Vos J.B., Rizzi A., Darrac D., Hirschel E.H. Navier —Stokes solvers in European aircraft design // Progress in Aerospace Sci. 2002. Vol. 38. P. 601-697.

[5] Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука. Сиб. отд-ние, 1967. 196 с.

[6] Белоцерковский О.М. Численное моделирование в механике сплошных сред. М.: Наука, 1984. 520 с.

[7] Годунов С.К., Забродин А.В., Иванов М.Я. и др. Численное решение многомерных задач газовой динамики. М.: Наука, 1976. 400 с.

[8] Самарский А.А. Теория разностных схем. М.: Наука, 1989. 614 с.

[9] 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. Р. 87-108.

[10] Briley W.R., McDonald H. Solution of the 3D compressible Navier — Stokes equations by an implicit technique // Lect. Notes in Phys. 1975. Vol. 35. Р. 150-178.

[11] Пейре Р., Тейлор Т.Д. Вычислительные методы в задачах механики жидкости. Л.: Гидрометеоиздат, 1986. 320 с.

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

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

[14] Ковеня В.М. Разностные методы решения многомерных задач. Курс лекций. Новосибирск: Новосибирский гос. ун-т, 2004. 146 с.

[15] Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1978. 736 с.

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

[17] Ковеня В.М., Слюняев А.Ю. Модификации алгоритмов расщепления для решения уравнений Эйлера и Навье —Стокса // Вычисл. технологии. 2007. Т. 12, № 3. С. 71-86.

[18] Ковеня В.М., Слюняев А.Ю. Алгоритмы расщепления при решении уравнений Навье— Стокса // Журнал вычисл. математики и математ. физики. 2009. Т. 49, № 4. С. 700-714

[19] Ксуемуа V.M. Splitting method in the problems of CFD // Comput. Fluid Dynamics J. 1999. Vol. 8, No. 2. P. 186-194

[20] Ковеня В.М. Алгоритмы расщепления при решении многомерных задач аэрогазодинамики. Новосибирск: Изд-во СО РАН, 2014. 280 с.

[21] Ковеня В.М., Тарнавский Г.А., Чёрный С.Г. Применение метода расщепления в задачах аэродинамики. Новосибирск: Наука. Сиб. отд-ние, 1990. 247 с.

Поступила в 'редакцию 22 января 2014 г., с доработки — 4 апреля 2014 г.