Модификация численного метода Годунова для уравнений газодинамики

Том ХЫН

УЧЕНЫЕ ЗАПИСКИ ЦАГИ

2012

№ 4

УДК 533.6

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

Ю. М. ЛИПНИЦКИЙ, А. В. САФРОНОВ

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

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

ВВЕДЕНИЕ

Условие неубывания энтропии наряду с выполнением интегральных законов сохранения является одним из необходимых условий применимости численных методов для расчета течений газа с разрывами. Решение энтропийной проблемы в численных расчетах впервые предложено Дж. фон Нейманом [1] путем введения искусственной вязкости, которая проявляется в градиентных зонах течения, «размазывая» разрывы.

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

Схематически начальный разрыв газа с различными состояниями в левом и правом полупространстве распадается на три волны: левую волну, контактный разрыв и правую волну. Обозначим их скорости И^, ф и Ид соответственно. Левые и правые волны могут быть в зависимости от начального перепада давления как волнами разрежения, так и волнами сжатия (скачками). Точное решение задачи Римана единственно, зависит от автомодельной переменной, скорости волн удовлетворяют неравенству ф < И < Заданной функции потока (компонентам вектора потока) отвечают два решения: одно соответствует дозвуковому случаю, другое — сверхзвуковому. Поэтому можно представить, что функция

ЛИПНИЦКИЙ Юрий Михайлович

доктор технических наук, профессор, начальник центра теплообмена и аэродинамики ЦНИИмаш

САФРОНОВ Александр Викторович

кандидат физико-математических наук,

начальник отдела газодинамики старта ЦНИИмаш

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

Точное решение задачи распада разрыва приводит к трудоемким вычислениям, кроме того, в ряде физических случаев теоретическое решение этой задачи затруднительно. В этой связи широкое распространение получили методы на основе приближенного решения задачи Римана. К данному классу можно отнести известные схемы LxF (Lax-Friedrichs) [3], Русанова [4], HLL (Harten, Lax, Leer) [5], Roe [6], EO (Engquist — Osher) [7], HLLC (Toro, Spruce, Speares) [8].

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

Решением проблемы энтропии в численном решении гиперболических уравнений является также применение кинетического метода релаксации, например [9 — 13]. В этом варианте исходная нелинейная гиперболическая система уравнений преобразуется в линейную систему с релаксационным источником. Причем кинетическая интерпретация численных методов упрощает их энтропийный анализ.

Сравнительный анализ численных методов на основе решения задачи Римана и других подходов, применяемых для уравнений газодинамики можно найти, например, в работах [14, 15].

Интенсивное развитие методов Годунова связано с повышением порядка вычислений на основе процедур интерполяции параметров к граням ячейки сетки с ограничителями, обеспечивающими монотонность схемы. Первый такой ограничитель предложен Колганом [16], изложение теории и обзор вариантов можно найти, например, в работах [17 — 19]. После интерполяции к граням ячейки (реконструкции) параметры на границе ячейки в схемах повышенного порядка вычисляются из задачи распада разрыва аналогично схеме 1-го порядка. Поэтому повышение порядка не устраняет принципиальные недостатки и в первую очередь необходимо выполнение условия неубывания энтропии в базовых численных схемах 1-го порядка, которые анализируются в настоящей работе.

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

Из рассматриваемых методов доказанным является выполнение условия неубывания энтропии в газодинамических расчетах по схемам LxF, Годунова, Русанова, HLL, EO, а также HLLC в кинетической интерпретации [12, 13] и в массовых переменных [20].

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

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

с выпуклой функцией потока/(и). Для определенности примем /"(и) > 0 .

Уравнение (1) может иметь разрывные решения (скачки), скорость распространения которых 5 определяется соотношением Ренкина — Гюгонио: /(и+) - /(и_) = s(u+-и_). В случае гладких решений функция и(х, ¿) постоянна вдоль характеристик Жх/Ж = /' (и).

Энтропийное условие допустимости разрыва Олейник [22] заключается в том, что, например, в случае 5 > 0, характеристическая скорость за разрывом должна быть больше скорости распространения разрыва, а характеристическая скорость перед разрывом — меньше:

1. ЭНТРОПИЙНЫЙ АНАЛИЗ РАЗНОСТНЫХ СХЕМ

Щ + fx = 0

u (x,0) = u0(x), 0 < x < 1, 0 < t < T,

f'(u- ) > 5 > f'(u+ ).

Запишем разностную схему для уравнения (1) в консервативном виде:

Mf+1 = un — с

At л

где с = — ; n — номер шага по времени с интервалом At; i — номер ячеики сетки по оси х c раз-Ax

биением Ax.

В методах типа Годунова поток на границе ячеек fi + 1/2 вычисляется из решения задачи распада разрыва с начальными данными, соответствующими параметрам в соседних ячейках сетки (i, i + 1). При этом схема устойчива при числе Куранта CFL = smaxc < 1, где smax — максимальная скорость распространения возмущений в расчетной области.

Как упомянуто, кроме выполнения законов сохранения для единственности численного решения необходимо выполнение энтропийного закона. Условие, обеспечивающее неубывание энтропии в ячейках сетки при численных расчетах, являющееся следствием теоремы Олейник, введено Ошером [21]. Это условие (E-flux), которое необходимо выполнять при аппроксимации потока на границе ячеек (fi+1j2) имеет следующий вид:

f+1/2 — f

иг+1 — щ

Поток в методе Годунова на границе ячеек (1/,+°д,унов ), который вычисляется из точного решения задачи распада разрыва с начальными данными, соответствующими параметрам в соседних ячейках сетки, может быть представлен в виде [21]:

Я

i+1/2

Годунов =

uc[ui ,ui+1]

min f (u), u < u+1,

uc[u. ,u i+1 ]

' i+1

max (f (u)),u > u,+1.

Поток в методе Годунова удовлетворяет условию Е-/их, причем

/Е * /Годунов, щ < Щ+1,

/Е * /Годунов, Щ > иг+1.

С целью сравнительного анализа диссипативных и энтропийных свойств численных методов разностное уравнение (3) представляется в следующем виде:

иГ1=ип _-/- /л+2 а+1/2

где Q —схемная вязкость

Й+1/2 =с

^ fi+1 + fi 2fi+1/2 ^

u,-,i — u,

Н+1

В этом случае аппроксимация потока на границе ячеек выражается следующим образом:

¿+1/2 = /Ш2+/- _ 2- е1+1/2

2 2-

Схемная вязкость в методе Годунова имеет вид [21]:

Q

Годунов _

+1/2

= с max

uc[uj ,Uj+1

( f+1 + f - 2f (u) > ui+1 - ui

Энтропийному условию удовлетворяют методы, в которых при аппроксимации потока на границе ячеек схемная вязкость больше схемной вязкости метода Годунова [11]:

П > пГодун°в Q+1/2 г q+1/2 ,

что аналогично выполнению условия E-flux (4).

Согласно [11] рассматриваемые численные схемы монотонны и устойчивы при выполнении условия CFL и в случае

0+1/2 < Q+1/2 * QLF/2 - 1.

Здесь =Сfi+1/2 — схемная вязкость в методе Roe [7], Пил/2 — схемная вязкость

в методе Лакса LxF [3], принимающая максимальное значение.

Приведем сравнительный анализ величин схемной вязкости методов на основе приближенного и точного решений задачи распада разрыва. Ниже рассматриваются параметры на границе ячейки сетки. Индекс «L» относится к i-й ячейке сетки, индекс <R» — к i + 1-й, индекс «i + 1/2» — к параметрам на границе этих ячеек.

Схемная вязкость в методе EO [7, 21] имеет вид:

uR

Л f '(u)| du

qeo Qi+1/2

= с

uL

u p u :

В работе [21] показано, что поток EO удовлетворяет условию E-flux.

«Двухволновые» схемы HLL, Русанова и LxF представляются из соотношений Ренкина — Гюгонио, записанных для двух скачков (рис. 1):

(f - fL ) = sL (u*~ uL X (f - fR ) = sR (u - uR ).

Здесь SL, — скорости левой и правой волн соответственно. Причем для рассматриваемой выпуклой функции потока < ¿к.

Аппроксимация потока на границе ячеек из соотношений (14) имеет вид:

f = f * = fLsR - fRsL + sLsR

fi+1/^ — J -

при этом

i+1/2

= u =

sR - sL

uRsR uLsL fR + fL

sd st

Для определения скоростей волн используются следующие варианты:

•л =-$ь = 1/о,

зь = шт

В случае (17) поток (15) представляет схему ЬхГ.

Рис. 1. Геометрия энтропийного E-flux анализа

В этой схеме используется максимальная оценка скоростей волн по всей расчетной области. Заметим, что в этом случае поток f (u*) из (16) отвечает схеме 2-го порядка LxW [23].

При вычислении скоростей волн согласно (18) поток (15) соответствует схеме Русанова, в которой применяется локальный максимум скоростей волн. При выборе варианта (19) аппроксимации потока на границе ячеек (14) соответствуют схеме HLL.

Рассмотрим случай смены знака производной функции f (u) при ur > ul. Этот случай является аналогом волны разрежения и характерен для исследования энтропийных проблем. Геометрическая интерпретация потоков рассматриваемых схем представлена на рис. 1.

В соответствии с (5), (6) поток на границе ячеек схемы Годунова соответствует минимуму функции f (u) на рассматриваемом интервале, а схемы, удовлетворяющие условию E-flux, находятся ниже этого уровня, как показано на рис. 1. В данном случае поток EO (13) совпадает с потоком Годунова. Поток HLL находится в зоне E-flux, на пересечении касательных к функции f (u), проходящих через точки (ul, fL) и (ur, JR). Потоки Русанова и LxF находятся также в зоне E-flux ниже потока HLL, поскольку линии, проведенные из точек (ul, fL) и (ur, JR) в этих схемах ближе к вертикали, чем в методе HLL. Очевидно, выбор скоростей волн (17) — (19) обеспечивает выполнение условия неубывания энтропии в двухволновых схемах LxF, Русанова и HLL, для любой выпуклой функции f(u).

В случае изменения знака характеристик (fL < 0, fR > 0), поток Roe, как показано на рис. 1,

не удовлетворяет энтропийному условию. На рис. 1 показан поток f (u*), которой при определении u по формуле (16) с выбором скоростей волн в виде (17) представляет собой схему LxW [23], при этом поток f (u*) больше потока схемы Годунова, но меньше потока Roe. В рассматриваемом случае чем меньше поток на границе ячеек, тем больше схемная вязкость (8), поэтому при смене знака характеристик схема LxW 2-го порядка предпочтительнее схемы Roe.

Согласно [11] и в соответствии с данными рис. 1 известные потоки по величине схемной вязкости располагаются следующим образом:

s^Roe ^ .^Годунов ^ глЕО ^ s^HLL ^ ^Русанов ^ ,^¡LxF _ 1 Qi+1/2 < Qi+1/2 < Qi+1/2 < Qi+1/2 < Qi+1/2 < Qi+1/2 = 1

На рис. 1 показан способ учета контактного разрыва в задаче распада разрыва с выполнением условия E-flux. Алгоритм заключается в следующем: по параметрам u , полученным с помощью формулы (16), находится наклон касательной в этой точке s* = f' (u*), затем вычисляется

поток f(u ), и далее определяется искомое решение, находящееся на пересечении данной касательной линии с линиями, обозначающими левую и правую волны. Потоки, отвечающие этому

обозначены f и Jr , а переменные ul и ur соответственно. Для этого случая можно записать следующие соотношения на разрывах:

При этом в силу выпуклости функции f(u) выполняются неравенства

Sl < s* = f'

Выбор аппроксимации на границе ячеек зависит от знака s* :

fi+1/2 = fL при s* > 0 и¿+1/2 = fR при s* < 0.

Решение задачи на основе соотношений (20) можно представить в виде:

fi+1/2 = f* [f

-Sr

У S* S l SR S * j

На рис. 1 данные потоки обозначены XContact, X — означает способ выбора скоростей волн sl и Sr, соответствующий схемам LxF, Русанова, HLL (17) — (19). Из рис. 1 видно, что одно из

пересечений касательной к функции f (u) в точке (u*, f (u*)) с линиями, обозначающими левую и правую волны, при выпуклой функции f (u) в случае f < 0, fR > 0 всегда будет лежать ниже минимума этой функции. Поэтому потоки XContact удовлетворяют энтропийному условию E-flux.

Из анализа соотношений (22),(23) и данных рис. 1 выявлена следующая возможность модернизации двухволновых схем. Если в качестве скорости средней волны s* принять предельные

значения, т. е. из точки (u*, f(u*)) провести линии, параллельные левой и правой волне, то получим следующее значение параметра ю в уравнении (22):

ю=ю

GFORCE

= min

-Sr

\ sr s

L JR

Sd - S

L

™=™contact= min

При выборе sr = -sl = smax, полученная схема (15), (16), (22), (24) соответствует схеме GFORCE [24], поэтому в соотношении (24) параметр ю обозначен &gforce.

Новые схемы (15), (16), (22), (24) с выбором скоростей волн (18) или (19) являются обобщением схемы GFORCE. На рис. 1 положение данных потоков обозначено XGFORCE. Из геометрии видно, что полученные схемы удовлетворяют условию E-flux. Потоки XGFORCE, основанные на максимальной оценке скорости средней волны (рис. 1), которая всегда выше скорости контактного разрыва, «размазывают» контактный разрыв, но их схемная вязкость (8) меньше, чем соответствующих потоков при ю = 0 (LxF, Русанов, HLL).

Исходя из вышеизложенного, предложенные модификации численных методов (XContact и XGFORCE) по величине схемной вязкости можно расположить следующим образом:

ÓRoe ^ ^Годунов ^ r>EO ^ ^XContact ^ ^XGFORCE ^ s^HLL ^ ^Русанов ^ s^LxF _ л i+1/2 < Qi+1/2 < Qi+1/2 < Qi+1/2 < Qi+1/2 < Qi+1/2 < Qi+1/2 < Qi+1/2 = 1

В соответствии с (12) схемы XContact и XGFORCE устойчивы и монотонны.

В общем случае при выборе параметра ю в диапазоне 0 <ю <ю contact схемы (15), (16), (22) с выбором скоростей (17) — (19) удовлетворяет условию E-flux, поэтому возможны и иные варианты, с различными переключателями.

2. УРАВНЕНИЯ ГАЗОДИНАМИКИ

Рассмотрим двумерный случай в расщеплении по оси х. Уравнения газодинамики представляются в следующем виде:

и + Ех = 0,

где и — вектор консервативных переменных и Е — вектор потока:

р pw

U = pw , F = pw2 + P

pv pwv

PE _ pEw + Pw

Здесь р — плотность газа; V, V — нормальная и касательная к грани ячейки компоненты скорости; Е = е + м> 2/2 + V 2/2 — полная энергия на единицу объема; е = р/р/(у - 1) — внутренняя энергия; р — давление; у — показатель адиабаты.

Методы, описанные для скалярного закона сохранения в разделе 1 (3), (15), (16), (22), (23), (24), обобщаются для системы (25), путем замены величины и на вектор консервативных переменных и, величины/— на вектор потока Р, а значения скоростей волн SL, я* и sя соответственно на Жь, Ж* и Жя:

тП+1

где

иг1 = ип _-

Жя _ Жь ияЖя _ иьЖь _ Ря + рг

и* = -я"я

Жя _ Жь

Р,+1/2 = Р [Р (и ) _ Р ],

_ЖЬ

Жя

Ж* _ Жь Жя _ Ж *

или

ю = юОРОяСЕ = т1П I

_ЖГ

Жя

Жя _ Жь Жя _ Жь ,

где

Поток Р (и *) в формуле (28) вычисляется следующим образом:

Р (и *) = Ж,и * + 0*,

е* = {0,/,0,/Ж*}7 .

Скорость (Ж*) и давление (р*) на контактном разрыве определяются из соотношений на разрывах для уравнений газодинамики в лагранжевых переменных [20] с введением массовых скоростей через поверхности левой ть = рь(^ь - Жь) и правой тя = ря(Жя - ^я) волн:

Ж* = = wя =

* шя^я + т1Уь _Ря + Рь

ть + Шя

Р = Рь = Ря =

тьРя + тяРь _ тьтя

т

т

Уравнения в универсальном виде (28), полученные из соотношений на разрывах, в зависимости от выбора скоростей волн Жя, Жь и способа вычисления параметра ю включают в себя девять численных методов: классические схемы ьхР, Ньь, Русанова (ю = 0), их варианты с учетом контактного разрыва Х-СоШаМ (29) и модификаций Х-ОРОяСЕ (30), где «X» обозначает соответствующий способ выбора скоростей волн.

Выбор скоростей волн:

1) Жя = _Жь = 1/ - , соответствует схеме ьхР;

2) Жя = _Жь = max(abs(wя) + ся, abs(wL) + сь), соответствует схеме Русанова;

ю = Юсоп(ас1 = т1П

3) Ж = тт(м>ь - е^, - ер), Жр = шах(^ + е^,, ^^ + ер), соответствует схеме ИЬЬ.

Следует подчеркнуть, что все перечисленные подходы энтропийно обоснованы.

Заметим, что при аппроксимации потока на границе ячеек сетки по формулам (26) — (33) не требуется привлечения уравнения состояния. Это упрощает применение подходов для сложного уравнения состояния [25, 26] и расчетов многокомпонентных течений с переменными теплофизическими свойствами [27].

2 4 6 8 10 12 14 16 18

Рис. 2. Картина заполненных контуров плотности (результаты расчетов по схеме ЯЬЬ-СоШае^

2 1.6 1.2 0.8 0.4

Рре

х /Ra

rms -1 -

-2

-3

-4

-5

-6

12

16

20

RK3, KorerS

----HLL-Cortact

Годунов —  — HLL-GFORCE

Л .

Рис. 3. Изменения давления на оси струи

III liter

4000 8000 12000 16000 20000

Рис. 4. Сходимость расчетов

Обобщение изложенных методов на трехмерный пространственный случай аналогично методу Годунова [2].

В качестве иллюстраций на рис. 2 — 4 даны результаты численных расчетов сверхзвуковой осесимметричной струи воздуха, истекающей в затопленное пространство. Уравнения Эйлера решались методом установления. Параметры на срезе сопла: число Маха Ма = 2, отношение давления на срезе сопла к внешнему давлению ра/ре = 1.5, температура торможения Т00 = 300 К. Струя истекает справа налево, срез сопла расположен от нуля до 1 по вертикальной оси. На срезе сопла задавался сверхзвуковой поток, на оси струи — условия симметрии, на гранях ячеек, примыкающих к свободным границам, потоки задавались из решения задачи Римана с параметрами, соответствующими состояниям в примыкающей ячейке и невозмущенной внешней среде. Расчетная область занимала 3 х 20 радиусов выходного сечения сопла (Ра = 1). Использовалась квадратная сетка с количеством разбиений на радиус среза сопла Та = 20.

Расчеты проводились по ТУБ схеме Рунге — Кутта 3-го порядка по времени с ограничителем Корена 3-го порядка по пространству (РК3-Когеп) [19] с решением задачи распада разрыва по схемам Годунова: с точным решением задачи Римана, а также представленными ЯЬЬ-СоМаМ и ЯЬЬ-ОЕОРСЕ.

В потоке образуются скачки, волны разрежения и контактный разрыв (рис. 2). Разрывы проходят поперек ячеек сетки и рассчитываются сквозным образом (без выделения).

0

4

8

0

На рис. 3 приведены изменения давления на оси струи, обозначения кривых соответствуют рис. 4.

На рис. 4 показана сходимость по времени в виде изменения логарифма максимального (в расчетной области) приращения плотности за шаг счета в зависимости от количества шагов счета по времени (итераций).

Из рис. 3, 4 видно, что результаты расчетов с применением предложенных способов решения задачи распада разрыва близки к схеме Годунова.

ЗАКЛЮЧЕНИЕ

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

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

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

Результаты работы могут быть применены для решения широкого круга задач аэрогазодинамики.

Работа выполнена при поддержке грантов РФФИ № 10-08-00501, № 10-01-00677.

ЛИТЕРАТУРА

1. Neumann J., Richtmyer R. A method for the numerical calculation of hydrody-namic shocks // J. Appl. Phys. 1950, 21, № 3, р. 232 — 237.

2. Численное решение многомерных задач газовой динамики / Под ред. С. К. Годунова. — М.: Наука, 1976.

3. Lax P. D. Weak solutions of nonlinear hyperbolic equations and their numerical computation // Comm. Pure Appl. Math. 7, 159 (1954).

4. Русанов В. В. Расчет взаимодействия нестационарных ударных волн с препятствиями // ЖВМ и МФ.1961, 1, N 2, р. 267 — 279.

5. Harten A., Lax P. D., B. Van Leer. On upstream diffrencing and Godunov-type schemes for hyperbolic conservation laws // SIAM Review. 1981. 25, N 1, р. 35 — 61.

6. R o e P. L. Approximate Riemann solvers, parameter vectors, and difference schemes // J. Comput. Phis. 1981. 43, N 2, р. 357 — 372.

7. Enguist B.,Osher S. One-sided difference approximation for nonlinear conservation laws // Math. Comp. 36(1981), р. 321 — 351.

8. Toro E. F., Spruce M., Speares S. Restoration of the contact surface in the HLL Riemann solver // Shock Saves. 1994, 4, р. 25 — 34.

9. Jin S., X i n Z. P. The relaxation schemes for systems of conservation laws in arbitrary spacedimensions // Comm. Pure Appl. Math. 1995. 48 : 235 — 276.

10. L e V e q u e R. J., P e l a n t i M. A class of approximate riemann solvers and their relation to relaxation schemes // J. Comput. Phys. 2001, 172, р. 573 — 591.

11. Tadmor E. Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems // Acta Numerica (2003), p. 451 — 512. Cambridge University Press, 2003.

12. Bouchut F. Entropy satisfying flux vector splittings and kinetic BGK models// Nu-mer. Math. 2003, 94, р. 623 — 672.

13. Сафронов А. В. Кинетические схемы для уравнений газодинамики // Вычислительные методы и программирование. 2009. Т. 10, с. 62 — 74.

14. Toro E. F. Riemann solvers and numerical methods for fluid dynamics. — Berlin. Heidelberg: Springer-Verlag, 2009.

15. Воронич И. В. Сравнительный анализ группы численных методов газовой динамики. — М.: МФТИ, 2007, 151 с.

16. Колган В. П. Применение принципа минимальных производных к построению конечно-разностных схем для расчета разрывных решений газовой динамики // Ученые записки ЦАГИ. 1972, Т. III, № 6, с. 68 — 77.

17. Bram van Leer. Upwind and high-resolution methods for compressible flow: From Donor Cell to Residual-Distribution Schemes // Commun. Comput. Phys. 2006. V. 1, N 2, p. 192 — 206.

18. B e r g e r M., A ft o s m i s M. J. Analysis of slope limiters on irregular grids // AIAA Paper 2005-0490. 2005.

19. Сафронов А. В. Оценка точности и сравнительный анализ разностных схем сквозного счета повышенного порядка // Вычислительные методы и программирование. 2010, 11.

20. Сафронов А. В. Разностный метод для уравнений газодинамики из соотношений на разрывах // Математическое моделирование. 2008. 20, № 2, с. 76 — 84.

21. O s h e r S. Riemann solvers, the entropy condition, and difference approximations // SIAM J. Numer. Anal. 1984; 21(2) : 217 — 235.

22. Олейник О. А. О единственности и устойчивости обобщенного решения задачи Коши для квазилинейного уравнения // Успехи мат. наук. 1959. 14, № 2 (86), с. 165 — 170.

23. Lax P., Wendro ff B. Systems of conservation laws // Comm. Pure Appl. Math. 13 (1960), р. 217 — 237.

24. Toro E. F., Titarev V. A. MUSTA fluxes for systems of conservation laws // J. of Comp. Phys. 216 (2006), р. 403 — 429.

25. Сафронов А. В. Разностный метод решения уравнений гидродинамики с двучленным уравнением состояния на основе энтропийно-согласованного приближенного решения задачи Римана. Высокопроизводительные вычисления в задачах механики и физики: Сб. научных трудов, посвященный 75-летию со дня рождения А. В. Забродина. — М.: ИПМ им. Келдыша, 2009, с. 143 — 149.

26. Липницкий Ю. М. , Сафронов А. В. Численные методы Годунова типа для уравнений газодинамики (Риман — солверы). — Тезисы докладов // Третья всероссийская конференция «Вычислительный эксперимент в аэроакустике». — Светлогорск, Калининградской области. 2010.

27. Сафронов А. В. Численный метод расчета струй продуктов сгорания при старте ракет // Космонавтика и ракетостроение. 2007. № 1(46), с. 72 — 79.

Рукопись поступила 9/III2011 г.