Сравнение точности и сходимости некоторых TVD-схем Текст научной статьи по специальности «Математика»

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

Том 5, № 5, 2000

СРАВНЕНИЕ ТОЧНОСТИ И СХОДИМОСТИ НЕКОТОРЫХ TVD-СХЕМ*

Д. В. Чирков Новосибирский государственный университет, Россия

С. Г. ЧЕРНЫЙ Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: [email protected]

Versions and modifications of the upstream schemes suggested by Harten, Yee and Chakravarthy have been considered. An attempt to employ a unified approach to Chakra-varthy's principle of the TVD-schemes construction and to the systematization of the well-known limiter functions has been made. One of the main objectives of this study is to investigate the applicability of the TVD-algorithms to the solution of stationary problems using the relaxation method. It is shown that those TVD-schemes which have the highest accuracy of discontinuity representation cannot provide the residual convergence to the macine zero. The optimum choice combining satisfactory accuracy with the high rate of convergence is Chakravarthy's scheme with Van Leer's limiter.

Введение

Целью данной работы является сравнение различных TVD-схем на предмет их точности и скорости сходимости к решению стационарной задачи, отыскиваемому методом установления. Необходимость подобного исследования назрела в процессе разработки универсального численного алгоритма, пригодного для расчета течений газа в широком диапазоне скоростей: от малых дозвуковых до гиперзвуковых. Как известно, уменьшение характерного числа Маха течения приводит к заметному ухудшению сходимости методов установления, что связано с возрастающим при M ^ 0 различием между минимальным и максимальным по модулю собственными значениями матриц Якоби векторов потоков. В работах [1, 2] предложен достаточно эффективный метод ускорения сходимости — метод пред-обуславливания (preconditioning), заключающийся в модификации члена с производными по времени в областях, где локальное число Маха меньше единицы. Такая модификация может быть проведена для любого метода установления, а эффективность полученного универсального алгоритма будет в первую очередь зависеть от свойств базовой схемы. В частности, базовая схема должна обладать достаточной точностью и высокой скоростью сходимости к стационарному решению при расчете сверхзвуковых течений. Наиболее приспособленными для расчета таких течений и получившими в настоящее время широкое

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №98-01-00742, и Федеральной целевой программы "Интеграция", проект 274.

© Д. В. Чирков, С. Г. Черный, 2000.

распространение являются ТУБ-схемы. Поэтому поиск эффективной базовой схемы в данной работе будет проводиться в классе ТУБ-схем. Существует немало работ, посвященных сравнению точности различных ТУБ-схем [5, 7-10], в то время, как их сходимость очень редко служит объектом пристального исследования. Дальше всех в этом отношении продвинулись отечественные авторы Ю.М. Белецкий и др. (см. [7]), детально сравнившие варианты явных схем Хартена и Чакравати при расчете различных режимов течения в сопле Лаваля.

В предложенной работе мы ограничимся рассмотрением неявных противопоточных схем типа Хартена [3-6, 9] и типа Чакравати [11-13], имеющих первый порядок аппроксимации по времени и второй по пространству, оставив за рамками широкий класс МиБСЬ-схем, также представляющих определенный интерес.

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

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

1. Построение ТУВ-схем

1.1. Явные схемы для скалярного закона сохранения

Для скалярного уравнения

ди + д/{и) = 0 дЬ дх

рассмотрим явную разностную схему

АЬ

= и» - А/+1/2 - /7-1/2), А =—, (2)

где /±1/2 — потоки, вычисленные на слое п. В дальнейшем верхний индекс п, обозначающий величины, взятые на п-м слое, опускается. В общем случае /+1/2 для схемы первого порядка с разностями против потока имеет вид

А7+1/2« = «7+1 - «7,

/7+1/2 = т/+1 + /7 - ^{а7+1/2)А7+1/2«], (3)

ч, N > е,

Ф{*)={ г2 + е2 ,, (4)

-, Ы< е.

2е

Слагаемое ^{а7+1/2)А7+1/2м отвечает за численную вязкость схемы первого порядка. Обращение в нуль этого члена может приводить к нефизическому решению, в котором нарушается энтропийное неравенство, в частности, к появлению скачка разрежения. Для борьбы

1

с такого рода явлениями введена функция ф^), обеспечивающая энтропийную коррекцию (см. [3]).

Существует два принципиально различных подхода к построению противопоточных ТУБ-схем повышенного порядка аппроксимации по пространству. В первом из них, предложенном и детально проработанном Хартеном, в схему первого порядка с выражением (3) для потока вместо истинной потоковой функции f подставляется модифицированная функция f+д, где д выбирается таким образом, чтобы результирующий поток был второго порядка аппроксимации по пространству по отношению к функции f:

f3+1/2 = ^из+г + fj + д,-+1 + 9з - ф(а3+1/2 + Т.+1/2)Д3+1/2и],

{(Л+1 - f3 )/А3+1/2и, Д3+1/2и = О, (!) .■ « = 0,

_ /(д3+1 - д.)/Дз+1/2и, Д.+1/2и = 0, 7з+1/2 = \0, Дз+,/2« = 0. (5)

Для случая противопоточной схемы второго порядка по пространству и первого порядка аппроксимации по времени д. и д3+1 определяются следующим образом:

gj = minmod gj+1 = minmod

где

2 ^(aj+1/2)Aj+i/2u, 1 ^(aj_1/2 )Aj_1/2U

2 ф(aj+3/2)Aj+3/2u, 1 Ф(aj+1/2)Aj+1/2U

minmod(x,y) = sgn(x) ■ max(0,min(|x|,sgn(x)y)). На практике чаще используют альтернативный способ вычисления gj и gj+1:

gj = 2 ^(aj+1/2)gj, gj+1 = 1 Ф(aj+1/2) 9j+1, gj = minmod[Aj_1/2, Aj+1/2], gj+1 = minmod[Aj+3/2, Aj+1/2].

В этом случае поток fj+1/2 приобретает следующий вид:

f7+1/2

2

fj+1 + fj + 2 ^(aj+1/2)(gj+1 + gj) - Ф (aj+1/2 + Yj+1/2)Aj+1/2U

(8)

Исчерпывающее описание этого подхода, а также другие варианты вычисления gj можно найти в работах [3 - 5].

Рассмотрим подробно второй (почти классический) способ конструирования противопоточных TVD-схем повышенного порядка, который был использован в работах Чакрава-ти — Ошера [12], Роя [11] и многих других авторов. Назовем его подходом Чакравати.

Для начала рассмотрим случай, когда в (1) f = au, а = const, а > 0. Тогда поток, отвечающий противопоточной схеме второго порядка аппроксимации по пространству и первого по времени можно записать в виде

1 а fj+1/2 = 2 [fj+1 + fj - 1 a|j+1/2Aj+1/2U] + 2Aj_1/2U. (9)

6

7

Здесь первая часть соответствует потоку первого порядка, а последнее слагаемое повышает порядок аппроксимации по пространству до второго и обычно называется антидиффузионным членом. Анализ устойчивости методом Фурье показывает, что явная схема (2) с потоками (9) абсолютно неустойчива.

ТУО-схема получается из (2) и (9) путем домножения антидиффузионного члена в (9) на соответствующую функцию-ограничитель:

1 а

Л7+1/2 = 2 [Л'+1 + Л - N7+1/2^+1/2«] + 2А^-1/2иф-. (10)

Следуя методике Ши и Торо [14], выведем достаточные ограничения на функцию ф7-, при которых схема (2) с потоками (10) будет ТУО.

Воспользуемся следующим достаточным признаком: если разностная схема представи-ма в виде

«П+1 = Ц? - ^7-1/2 А7-1/2И + £7+1^7+1/2«

и

Я?_1/2 > 0, £+1/2 > 0, ^+1/2 + £+1/2 < 1, (11)

то она является схемой ТУО [3]. Введем обозначения:

¿7 = Д^, = (12)

Подстановка (10) в (2) дает

иТ1 = Ц7 - Аа ( Ц7 + ^з_1/2ифз

Ц7 - Аа + 1д7_1/2цф7 - и_1 - 1д7_з/2цф7_1^ = Ц7 - Аа 1 + 1 ф7 - 1Д7_1/2

Для удовлетворения (11) достаточно выполнения

0 < Аа ( 1 + 1 ф,- - 2< 1,

что равносильно

-2 < (ф. - ¿^ф^) < ^^, (13)

где с = Аа — число Куранта. Пусть ф7- > 0, тогда достаточным условием выполнения (13) будет

< 2(1 - с)/с,

Ф7 < 2^7, (14)

>7 > 0.

В большинстве предложенных ТУО-схем функция-ограничитель ф7- зависит лишь от ¿7:

Фз = Фз (¿7 )• (15)

Причем логично положить ф7- = 0 при ¿7 < 0, что соответствует локальному экстремуму в точке ]. Неравенства (14) определяют ТУО-область, в которой может находиться график

Рис. 1. Функции-ограничители. Верхняя граница ХУБ-области соответствует с = 0.4.

функции фз (6З). Заметим, что граница ХУБ-области зависит от с. Необходимое условие фз (1) = 1 следует из соображений сохранения второго порядка схемы (2), (10) на гладких решениях.

Ниже приведены наиболее известные функции-ограничители (см. также рис. 1). Ограничитель ЫШИСБ

фз =штшо^1,07), (16)

где функция minmod определена в (6). Ограничитель SUPERBEE

фз = шах[0, шт(1, 26з), шт(6З, 2)]. Ограничитель Ван Лира [15]

фз =

63 + |6з |

(17)

(18)

(19)

1 + 9з

Ограничитель ЫШИСБ с дополнительным сжатием

фз = [1 + ывз] • шinшod(1, 6з),

вз = З2и - А;-1/2ц|,, ы > 0.

|Дз + 1/2и| + |аз-1/2 и|

Очевидно, и в этом случае фз зависит лишь от 6З, т. е. всю комбинацию в правой части (19) можно рассматривать как ограничитель. В таком подходе нет необходимости оправдывать навешивание дополнительного сжатия на ЫШЫСБ, а затем доказывать, что оно не портит монотонности и точности схемы, как это делается, например, в [3, 4]. Легко убедиться, что при ы = 1 (19) совпадает с ограничителем Ван Лира. На рис. 1 приведен график функции (19) при ы = 2. Заметим, что если ы > 1, то этот график не лежит целиком в ХУБ-области.

Ограничитель Чакравати — Ошера

1_0 1 + g g_0

ф = —^minmod(1, ) + —^minmod(0/, в), 1 < в < ц—g• (20)

Здесь параметр 0 (не путать с 0/ и 0/) определяет тип схемы и порядок ее аппроксимации. Например, 0 = —1 отвечает противопоточной схеме второго порядка. Схема (2), (10) с ограничителем (20) при 0 = 1/3 соответствует противопоточной TVD-схеме третьего порядка аппроксимации, предложенной Чакравати и Ошером в [12].

Использование конкретной функции ф/ вносит свои ограничения на величину допустимого числа Куранта с. Условие ф/ (1) = 1 подразумевает, что заведомо с < 2/3.

На практике вместо функции-ограничителя обычно используют так называемую ограниченную вариацию переменной u. Условие (15), а также линейность функции-ограничителя относительно операции умножения всех его аргументов на скаляр позволяют представить произведение ф/Aj-1/2n, фигурирующее в (10), в эквивалентном виде

фAj_i/2u = limiter (Aj-1/2u, Aj+i^u). (21)

Например, если ф/ = minmod(1,0/), то

ф/Aj-i/2u = minmod(Aj-i/2U, Aj+i^u).

Основным достоинством такого представления является то, что выражение, стоящее в правой части (21), определено для любых значений Aj-i/2u и Aj+i/2n, в то время, как 0/, необходимое для вычисления ф/, не определено при Aj-i/2u = 0.

Заметим, что если f = au, a = const < 0, то TVD-поток f/+i/2 в (2) для противопоточной схемы следует вычислять по такой формуле:

1 a fj+i/2 = 2 [fj+i + fj — Hj+i/2Aj+i/2u] — 2Ai+3/2u0j+i. (22)

Рассуждения, аналогичные случаю a > 0, приводят к следующим ограничениям на функцию ф/.+р

'$■+! < 2(1 — с)/с,

ф/.+i < 20/+i, (23)

j > 0.

Как и прежде, будем предполагать, что Ф/+ = 0/+i(0/+i).

1.2. Противопоточные TVD-схемы для нелинейного уравнения

В случае, когда f (u) — произвольная гладкая функция, a(u) = -f представляется в виде

du

a(u) = a+(u) + a-(u),

+ |a(u), если a(u) > 0, - I 0, если a(u) > 0,

a+ = a- =

I 0, если a(u) < 0, I a(u), если a(u) < 0,

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

f j+i/2 = 2[fj+i + fj — ^(aj+i/2)Aj+i/2u] + 2[a+-i/2AJ-i/2u — a—+3/2 A/+3/2u], (24)

где функция ф определена в (4),

1 [тг1-^, Aj+1/2« =

Aj+1/2M = 0.

2

Тогда TVD-поток /?+1/2 через грань j + 1/2 приобретает вид fj+1/2 = 1[/j+1 + fj - ^(aj+1/2>Aj+1/2u] + 1[a+_i/2AJ-1/2« ■ ф - aT+3/2Aj+3/2u ■ 0j+l], (25)

где

фj = limiter(1, Qj>, 0j+1 = limiter(1, 0j+1>, +

Qj = —-7-, Qj+1 = ——7-• (26)

aj-1/2 Aj —1/2U aj+3/2Aj+3/2u

Выражение (25) равносильно следующему:

Л7+1/2 = 1[/7+1 + Л" - ф(а3+1/2)А3+1/2и] + +1[Нт^ег(а+_1/2Д3_1/2и, а++1/2Д7+1/2и) - Н^ет^+з^+з^«, 7+1/2 Д7+1/2и)]- (27)

Замечание. Можно получить альтернативное выражение для потока /7+1/2, если О7 и 07-+1 в (25) вычислять по прежним формулам (12). В этом случае

/7+1/2 = 1[/7+1 + /7 - ф(а7+1/2)Д7+1/2«] +

+1 [а+_ 1 /2limiter(Д7_ 1/2и, Д7+1/2и) - Тз^Н^ет^+З^, Д7+1/2и)]- (28)

Еще один вариант представления потока получается, если в (28) вместо а+_1/2 и а_+3/2

взять соответственно a++1/2 и a-+1/2. Тогда

2

2 V 2 1

f7+1/2 = 1[/j+1 + fj - ^(aj+1/2>Aj+1/2u] +

+ ô ô^(aj+1/2>[limiter(Aj_1/2M, Aj+1/2U> + limiter(Aj+3/2U, Aj+1/2U>] +

+ 1 aj+1/2[limiter(Aj_1/2u, Aj+1/2U> - limiter(Aj+3/2U, Aj+1/2U>]^ • (29)

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

Покажем связь схемы (2) с потоком (29) c наиболее распространенным вариантом схемы Хартена (2), (8). Если в качестве функции-ограничителя limiter используется minmod, то

limiter (A j_ 1/2 u, Aj+1/2U> = gj, limiter(Aj+3/2U, Aj+1/2U> = gj+1,

где <7 и <7+1 определены в (7), а если еще |а^+1/2| > е, то

2«7+1/2 [11т^ег(Д1?_1/2И, Д,+1/2и) - 11ш11вг(Д^+з/2М, Д^+1/2м)] =

= -8Щп(а7+1/2Ь7+1/2Д7+1/2и,

3

где 77+1/2 вычисляется в (5) через (7). Пусть теперь 107+1/21 > 2е. Тогда

/j+1/2 = 1 [/j+1 + - (|aj+1/2| + sign(aj+i/2)7j+1/2)Aj+1/2u] + 11 aj+1/21 (gj + gj+l)- (30)

Из определения функции minmod следует, что

1 < |Д 7+1/2 1 < |Д7+1/2|, ЙЩп((7) = вЩП^+О-

Тогда |77+1/2| < 21«7+1/21- Это позволяет в указанных предположениях записать следующие равенства:

1 «7+1/21 + 8Щп(а7+1/2)7.7+1/2 = 1 «7+1/2 + 77+1/2 1 = ^(«7+1/2 + Т7+1/2)-Окончательно получаем, что при |а7+1/21 > 3е

/7+1/2 = 2[Л+1 + /7 - ^(«7+1/2 + 77+1/2)Д7+1/2и + 1 ^(«7+1/2)((- + ^7+1)> (31)

что, очевидно, совпадает с потоком (8) для схемы Хартена. Итак, для нелинейного скалярного уравнения ТУБ-схема Чакравати с потоком (29) совпадает со схемой Хартена (2),(8), за исключением, разве что, некоторых окрестностей точек, в которых а(и) = 0.

1.3. Неявные противопоточные ХУО-схемы

Методика построения неявных ТУБ-схем подробно изложена в [4]. Представим здесь лишь основные моменты. Рассмотрим чисто неявную схему

,га+1 | \ ( (•«,+ ! (•«,+! N

«Г1 + /+1/2 - /-1/2) = j (32)

где

/3+1/2 = 1[/j++11 + /Г1 - N7+1/2Aj+1/2«n+1] +

+ 2[(a+-1/2)ra+1Aj-1/2«n+10n+1 - (a7+3/2)7+1Aj+3/2«7+10j7++11]. (33)

В первую очередь нас интересует, при каких условиях эта схема будет TVD. Пусть / = аи, а = const > 0, тогда

а

/П+1 = 7+1 I _ . Д , «7+1ф7+1

/3+1/2 = а«з + 2 ' Дз-1/2« -

Подстановка в (32) дает

«Г1 + Да (1 + 1- 1 j+131) Дз-1/2«га+1 = j (34)

где (7 из (12). В [4] доказано, что для того, чтобы схема (34) была ТУБ, достаточно выполнения одного условия

С7_+1/2 - (1 +1ФП+1 -1771) > о, (35)

которое удовлетворяется, в частности, если

^ - 7 (36)

ф > 0. V У

Таким образом, в неявном случае на функцию Ф7+1 накладывается меньше ограничений и ТУБ-область для нее шире.

Для практической реализации неявная схема (32) с потоком (33) линеаризуется методом Ньютона, причем члены, повышающие порядок аппроксимации по пространству, берутся на п-м слое. Линеаризованная схема записывается в А-форме:

(Д + а+_1/2А7_1/2 + а_+1/2А7+1/2^ ДиП+1 = -(/7+1/2 - /7-1/2), (37)

где АМП+1 = м7+1 - и7, а а+_1/2 и а+_1/2 вычисляются на п-м слое. К сожалению, не удается строго доказать, что линеаризованная схема (37) сохраняет свойство неувеличения полной вариации численного решения.

1.4. Обобщение на случай системы гиперболических уравнений

Рассмотренную выше теорию легко перенести на случай гиперболической системы уравнений с постоянными коэффициентами

д и д и

— + А — = 0, А = ИБЬ, Б = Diag(A1 , ...,Ат). (38)

д£ дх

Здесь И — вектор из т неизвестных, Ь — матрица, строки которой являются левыми собственными векторами матрицы А, И, — матрица, столбцы которой являются правыми собственными векторами матрицы А, А1, ... , Ат — собственные числа системы. Переход к характеристическим переменным W = ЬИ позволяет расщепить систему на т независимых скалярных уравнений

д W д W

иг + Баг = 0 (39>

для каждого из которых строится скалярная ТУБ-схема. Для нелинейной системы гиперболических уравнений

дИ + ^ (40)

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

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

dU д U

Ж + Aj+1/2 IX

— + A7+1/2 — = 0. (41)

Затем определяются локальные собственные числа А1+1/2, ... , Ат+1/2, локальные матрицы левых и правых собственных векторов Ь7+1/2 и К7+1/2, приращения характеристических переменных 07+1/2 = Д7+1/2W = Ь7+1/2(и7+1 — и7), и все сводится к рассмотренной выше ситуации. Технология обобщения ТУБ-схем типа Хартена на случай нелинейной гиперболической системы уравнений описана в [3, 5]. В качестве примера приведем схему Чакравати — Ошера второго порядка (0 = —1 в (20)) в том виде, в котором она была представлена в работе [12]:

ип+1 = ига — Д (рП+1/2 — 7/2)'

ГП+1/2 = 1[Г(и7) + Г(и7+1) — (ДГ++1/2 — ДГ_+1/2)] + 1[ДГ_+з/2 — ДГ+_1/2].

Здесь

AF±+1/2 = Rj+1/2D±+1/2L7+1/2A7+1/2U,

AF_+3/2 = R7+3/2D_+3/2L7+3/2 A7+3/2U,

AF+_1/2 = R7_1/2D"_1/2L7_1/2A7_1/2U,

D± = 1(D ±|D|),

+

м:

а ограниченные потоки ÄFF_+3/2 и AF +_1/2 вычисляются следующим образо

AF-+3/2 = R7+3/2 ■ minmod[a_+3/2,^^-+1/2], (42)

AF+_1/2 = R7_1/2 ■ minmod[^+_1/2 ,в^7++1/2],

где а± = D±LVAvU, v = j + 3/2, j + 1/2, j — 1/2. Операция minmod от векторных аргументов берется покомпонентно, ее результатом также является вектор.

Необходимо заметить, что, как и в случае одного нелинейного уравнения, здесь возникает произвол в выборе величин, подвергаемых ограничению. Это могут быть следующие приращения: AF = RDLAU, а = DLAU, а = LAU или собственно AU. Таким образом, множество TVD-алгоритмов можно разделить на dF-, а-, а- и dU-схемы. Наиболее естественным является вариант "а", отвечающий потоку (27) в скалярном случае, вариант "а" — это аналог схемы с потоком (28). Вследствие (42) рассмотренную выше схему Чакравати — Ошера следует отнести к разряду а-схем.

Дальнейшие упрощения связаны с точкой вычисления матриц R, D, L, используемых при конструировании Fn+1/2. Например, в [12] все матрицы, необходимые для потока Fn+1/2, вычисляются в точке j + 1/2. В этом случае схема для системы гиперболических уравнений представляет собой аналог схемы (2) с потоком (29) для скалярного уравнения.

В следующем разделе мы рассмотрим конкретные TVD-схемы на примере двумернной системы уравнений Эйлера.

2. Численные методы для системы уравнений Эйлера

Система уравнений Эйлера, описывающая двумерное плоское движение сжимаемого невязкого газа в декартовой системе координат, имеет вид

дЧ + + д его) = 0

д* дх ду

(43)

где

Я

/р\

рм

Vе/

/ ри \

2

Е

(

е

р^ рм^

\

рм2 + р рм^ ' р^2 + р

\(е + р)м/ \(е +

Здесь р — плотность, м,^ — компоненты скорости, р — давление, е — полная внутренняя энергия на единицу объема. Нас будет интересовать стационарное решение системы (43), полученное методом установления.

Неявная консервативная схема для системы уравнений (1)

Ч77+1 - Я:

Еп+1 _ Еп+1 е 7+1 _ е 7+1

+ Е ¿+1/2 Г,_1/2 + е7+1/2 е7_1/2 = о

А* Ах Ду

после линеаризации методом Ньютона записывается в операторной А-форме: А1 + ДХ А_+1/2 А,+1/2 + ДХ Аг+_1/2А*_1/2 + Ду А_+1/2 Д7+1/2 + Ду А++_1/2 Д7_1/2

(44)

ДЧ77+1:

1

п N га п \

- АХ(Е ¿+1/2 - Е ¿_ 1/2) - — (е7+1/2 - е7_1/2),

1

Ау'

(45)

где

ДЯ77+1 = Я77+1 - ЯП, Аг+1/2 = (^ЕдО^)*

Б,

= 1/2 Б,+ 1/2Ь,+1 /2,

¿+1/2 1,2 Л 3,4

+ 1/2 = , А , А , А Ы/2, Аг+1 /2 = ^ \+1/2 = М ± С

А±+1/2 = 2(А ± А))*+1/2, ф(А),+1/2 = И,+1/2^ (^(А1 ),ф(А2),ф(А3),ф(А4)),+1/2 Ь,+1/2 = И,+1/2ф(Б)г+1/2Ьг+1/2.

Функция определена в (4). Аналогично вычисляются А_+1/2 и А+_1/2 с той лишь

й (д о(д) у '

разницей, что А7+1/2 = I ———— , и соответствующие собственные числа на грани

V д Ч ) 7+1/2

2 + 1/2 равны А^^1/2 = V, А3^1/2 = V ± с. Здесь м, V — компоненты скорости, с — скорость звука на рассматриваемой грани ячейки. Все необходимые величины на грани определяются методом осреднения по Рою [11] через значения консервативных переменных Ч в прилегающих к грани ячейках.

Конкретная ТУБ-схема определяется способом вычисления явных потоков в правой части (45). Поток для схем типа Чакравати удобно представить в следующем виде:

1

ЕП+1/2 = 2[Гг+1 + - ф(А)г+1/2А,+ 1/2Ч] - Л¥,+ 1/2.

(46)

Здесь первая часть отвечает за первый порядок, а член Л содержит монотонизированные поправки, обеспечивающие второй или третий порядок аппроксимации по пространству. Схема Чакравати третьего порядка [12].

1 _ П 1 I 0

Л ¿+1/2 = —(М1 — N2) + —(М2 — N1), (47)

М1 = К;+з/2ттт°^аг"з/2,ва_+1/2], ^ = К;_1/2ттт°^а+1/2,ва+1/2], М2 = Ri+l/2minmod[aг_hl/2,вa_+з/2], Nl =

"¿±-1/2 = ± |С|^+1/2^+1/2(^¿+1 —

Третий порядок обеспечивается выбором 0 = 1/3. Сжимающий параметр в задается из промежутка

1 < в < ^ ■

причем, чем больше в, тем в меньшем числе узлов происходит переключение схемы с повышенного порядка аппроксимации на первый. Мы рассмотрим три варианта этого алгоритма: Б2В3 (в = 3), Б2В2 (в = 2), Б2В1 (в = 1).

Схема Чакравати второго порядка с ограничителями Ван Лира

Л¿+1/2 = 2(^+з/2 ■ ^¿+з/2 — ^-1/2 * ^гТ_1/2)' (48)

где компоненты векторов "¿__з/2 и ^¿_1/2 вычисляются по формулам

-е -е + | -е -£ |

~_е = "¿+з/2 ' "¿+1/2 + 1 "¿+з/2 ' "¿+1/21 (49)

"¿+з/2 = __е + -е , (49)

"¿+з/2 + "¿+1/2

"+е • "+е +1 "+е • "+е I

„ +е = "¿-1/2 "¿+1/2 + 1 "¿-1/2 "¿+1/21 , ^

"¿-1/2 = +е + +е . (50)

"¿-1/2 + "¿+1/2

Схемы Хартена. Схемы типа Хартена [4, 6, 9, 10] имеют поток вида

*?+1/2 = 2^+1 + ^ + R¿+1/2Ф¿+1/2], (51)

а компоненты вектора Ф¿+1/2 вычисляются следующим образом:

Ф+1/2 = 2 ^(Ле+1/2)(^е+1 + ^) — ^(Ае+1/2 + 7^1/2 К+1/2, (52)

«¿+1/2 = ^+1/2(^+1 — / = 1 и Ле ) / (^+1 — ^е)/ае+1/2' ае+1/2 = 0

7¿+1/2 = 2 ^¿+1/2)|о, ае+1/2 = 0.

Рассмотрим три выражения для вычисления д*, использующие ограничители МШМОБ, Ван Лира и БИРЕИВЕЕ:

= а*+1Д (53)

* = а_1/2 ■ а*+1/2 + К_1/2 ■ а!+1/2| ( )

а* + а* , (54)

д* = 5 ■ тах[0, min(2|a*+l/2|,Sa*_l/2), Л^О^Д 25а*_1/2)], (55)

5 = 8Щп(а*+1 /2) •

Получаем три варианта схемы Хартена, которые назовем соответственно НагМШ, НагУЬ, НагБВ.

Гибридная схема Хартена — Йи. Заметим, что в случае линейной системы (38) теория позволяет нам использовать различные ТУБ-схемы для решения отдельных уравнений системы (39). В [4, 5] аналогично предложено поступать и в нелинейном случае. Авторы этих работ рекомендуют применять менее сжимающие ограничители (МШМОБ, Ван Лира) для тех компонент векторов д,, которым соответствуют собственные числа м + с и м - с, и более сжимающие, такие как, например, БИРЕИВЕЕ, для тех компонент векторов д,, которым соответствует собственное число м.

Примером такой гибридизации является противопоточная схема Хартена — Йи (см., например, [5]), на которую мы будем ссылаться в дальнейшем как на ИРУев. Здесь, если А1^1/2 = м, А3+1/2 = м ± с, то д,1'2 вычисляется по формуле (55), а д3'4 — по формуле (54).

Оригинальная схема Хартена. В первых работах Хартена по теории ТУБ-схем (см., например, [3]) был предложен следующий способ вычисления компонент вектора Ф,+1/2 в

(51): 1

где

= [1 + ^0M, д* = ^^о^^и^Н^^ ф(А*+1/2 )a¿+l/2], 0* = |а,+1/2 - а,_1/2| ш> 0 * |а*+1/2| + |а*_1/21 ' " '

Здесь, в отличие от схемы НагМШ, введено дополнительное сжатие, регулируемое параметром и minmod берется от величин типа а. Мы рассмотрим два варианта этой схемы: ОНагШ1 (^ = 1) и ОНагШ2 (^ = 2).

Симметричная схема. В данной работе для сравнения с предложенными противопо-точными схемами интересно будет рассмотреть одну из симметричных (или центрально-разностных, как принято в русскоязычной литературе) ТУБ-схем, детально исследованных Йи (см. [5, 6]). В этой схеме (назовем ее БУМ) поток 3$7+1/2 записывается в виде (51),

но

Ф*+1/2 = -ф(А*+1/2)(«*+1/2 - 5*+1/2) , (57) <5*+1/2 = 5 ■ max(0, т^^О^Д 25а*_1/2, 25а*+3/2, 25(а*_1/2 + а*+1 /2)) , (58)

5 = в^О^Д

Из представленных нами алгоритмов Б2Б1, Б2Б2, Б2Б3, Б2УЬ, 0НагШ1, 0НагШ2 относятся к а-схемам, а НагМШ, НагУЬ, НагББ, иРУее, БУМ — к а-схемам. Все они имеют первый порядок аппроксимации по времени и второй по пространству, за исключением Б2Б3 и Б2Б2, имеющих формально третий порядок аппроксимации по пространству.

В заключение скажем несколько слов о реализации алгоритмов. Стабилизирующий оператор в левой части (45) приближенно факторизуется на два оператора, один из которых содержит все разности, направленные назад, а другой — все разности, направленные вперед [16]:

(] - д-А++_1/2Тг_1 - — А+Т-_1) ■ ]х (59)

ДхА+_1/2Те_1 Ду^1/2^"^ ' 1 . 1

где

х (] + ДХ А"+1/2т^+1 + Ду A7+l/2Tj+l)ДQn+1 = > (60)

] = Д1 + ДХ А+_1/2 - ДХ А_+1/2 + ДУ А+_1/2 - ДУ А_+1/2' (61)

а Тт±1 — операторы сдвига на единицу вперед или назад по направлению т. Обращение факторизованного оператора осуществляется бегущим счетом за два шага.

Однако специфика тестовой задачи, выбранной для сравнения алгоритмов, позволяет обратить неявный оператор в (45) без приближенной факторизации. Действительно, если х — сверхзвуковое направление, то А_+1/2 тождественно равна нулю во всей расчетной области. Поэтому ДЧ™+1 можно вычислить за один проход области, используя процедуру векторной прогонки при каждом фиксированном г. Этот способ заведомо избавляет нас от возможного влияния приближенной факторизации на процесс сходимости алгоритма. Численные расчеты по схемам Чакравати показали, однако, что предложенная ЬИ-факторизация не искажает картину сходимости алгоритмов, поэтому расчет по схемам Б2Б1, Б2Б2, Б2Б3, Б2УЬ проводился с использованием приближения (59), в то время как для всех остальных алгоритмов была реализована процедура прямого обращения неявного оператора векторными прогонками.

Для реализации граничных условий за каждой границей вводятся два ряда фиктивных ячеек, значения искомых переменных в которых при £ = £га либо полагаются равными соответствующим параметрам набегающего потока, либо экстраполируются из внутренней части расчетной области, в зависимости от типа граничных условий. При переходе на следующий временной слой приращение Qnг'+1 — Q¿гj = ДQГj1 в фиктивных ячейках полагается равным нулю. Такой подход позволяет использовать общие процедуры вычисления потоков для всех внутренних ячеек расчетной области.

3. Результаты расчетов 3.1. Постановка тестовых задач

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

Тестовая задача 1 — это задача о прохождении косого скачка уплотнения через область. Можно рассматривать ее как фрагмент задачи об отражении скачка от твердой стенки. Расчетная область этого течения представляет собой квадрат [0,1] х [0,1] (ограниченный пунктирной линией на рис. 2), нижняя граница которого — твердая стенка. На левой входной границе держатся величины р\ = 1, щ = 1, У\ = 0, р = 0.084933, соответствующие состоянию перед фронтом ударной волны, а на верхней — р2 = 1.776135, и2 = 0.890755, у2 = -0.189218, р2 = 0.194178, характеризующие состояние за фронтом. Соответствующее число Маха на левой границе = 2.9, ^ = 30°.

® ^^^^

0 - 0)

у

□.□ 1.0 2.0 3.0 4.0

Рис. 2. Расчетная область задачи об отражении косого скачка уплотнения.

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

Во всех приведенных ниже расчетах шаг по времени Д1 полагался равным 0.1 (соответствующее число Куранта СЕЬ ~ 2), что близко к оптимальному значению временного шага.

3.2. Точность

На рис. 3, а, б, в представлены графики распределения плотности вдоль оси х (при у = 0.675), полученные для задачи 1 различными ТУБ-алгоритмами на сетке 20 х 20. Видно, что выделяются классы схем по точности:

(ИагЫШ, Б2Б1, ЯУЫ}, (иРУее, ОИагШ1, ИагУЬ, Б2УЬ, Б2Б2}, (ОИагШ2, Б2Б3, ИагЯБ}.

Для удобства восприятия результатов на каждом рисунке изображены графики плотности для схем отдельно взятого класса в сравнении с графиками плотности характерных представителей двух других классов. Из рисунков видно, что схемы второго и особенно третьего класса обладают вполне приемлемой точностью, размазывая разрыв на 6-7 узлов сетки. При этом следует иметь в виду, что фронт скачка наклонен под углом 30° к сечению, что заметно снижает точность воспроизведения разрыва по сравнению со случаем перпендикулярного к сечению расположения скачка. Схема OHarW2 дает слабый всплеск за скачком, возможно, это связано с тем, что ограничитель (19) не лежит полностью в ТУБ-области.

Рисунок 4 дает распределение скорости для некоторых характерных ТУБ-схем. Очевиден провал в профилях скорости сразу за разрывом (на удалении от скачка графики

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

3.3. Сходимость

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

Errn = max{ max {| RHS™ [k] |}},

1<fc<4

где — правая часть (45), вычисленная на n-м слое.

На рис. 5 и 6 представлены графики зависимости этой величины от n, характеризующие сходимость различных TVD-схем для задачи 1, полученные на сетках [20 х 20] и [15 х 20]. Во всех алгоритмах использовалась энтропийная коррекция вида (4), причем для схем D2B1, D2B2, D2B3 е = 0.01, а для всех остальных схем полагалось е = 0.05, если не указано противное.

Видно, что схема Чакравати при в =1 и схема HarMIN, совпадающие по точности, имеют и одинаковую скорость сходимости. Схожесть этих алгоритмов была замечена нами еще на этапе построения. Расчет по этим схемам требует около 300 итераций для достижения машинной точности. Симметричная схема SYM очень быстро достигает сходимости до значений Err ~ 10-10, которая в дальнейшем не улучшается.

Наиболее точные схемы D2B3, HarSB и OHarW2 обладают в целом неудовлетворительной сходимостью, причем глубина установления решающим образом зависит от выбранной сетки (сравните, например, OHarW2 на рис. 5 и OHarW2 на рис. 6). Как правило, при расчете по этим алгоритмам эволюция решения переходит в фазу незатухающих или даже периодических колебаний. Следует, однако, отметить, что для практических нужд сходимости схемы до значений Err ~ 10-4 оказывается вполне достаточно.

Схемы UPYee, HarVL, OHarWl занимают промежуточное положение как по точности, так и по сходимости. Необходимо отметить, что сходимость алгоритмов UPYee, OHarWl также существенно зависит от сетки: при расчете на сетке 15 х 20 решение, получаемое по этим схемам, выходит на периодические колебания с минимальной невязкой Err порядка 10-6.

3.4. Влияние энтропийной коррекции

Для большинства схем и входных параметров увеличение е в (4) наряду с ухудшением точности передачи разрыва влечет ускорение сходимости. Наиболее чувствительной к выбору параметра энтропийной коррекции оказалась схема Чакравати D2B2 на сетке 20 х 20 (рис. 7 и 8). Однако в случае, например, расчетов задачи 1 на сетке 15 х 20 по схемам D2B3, D2B2, OHarW1 увеличение е не позволяет вывести график истории сходимости из фазы периодических колебаний около значения Err = 10-4.

Рис. 3. Распределения плотности для задачи 1 в сечении у = 0.675.

Рис. 4. Распределения скорости для задачи 1 в сечении у = 0.675.

Рис. 5. Сходимость ТУБ-схем на задаче 1 (сетка 20 х 20).

Рис. 6. Сходимость ТУБ-схем на задаче 1 (сетка 15 х 20).

3.5. Полная задача об отражении скачка уплотнения

Рисунок 9 иллюстрирует распределения скорости в сечении у = 0.575, полученные для задачи 2 некоторыми характерными ТУБ-схемами на сетке 80 х 20. Очевидно, иерархия схем по точности сохраняется и при воспроизведении отраженного скачка. Обозначается область локализации нефизического "провала" в распределениях скорости (а также плотности), формирующегося за падающим скачком.

При этом глубина "провала" для алгоритма OHarW2 примерно вдвое меньше, чем для всех остальных схем. Из рис. 10, на котором представлена типичная картина изолиний

Рис. 7. Влияние параметра е в (4) на точность Рис. 8. Влияние параметра е в (4) на сходи-схемы Б2Б2 (у = 0.675, сетка 20 х 20). мость схемы Б2Б2 (сетка 20 х 20).

Рис. 9. Распределение скорости.

Рис. 10. Изолинии плотности для задачи 2, полученные по схеме Б2Б3.

плотности, видно, что "провал" представляет из себя вторичную волну, берущую начало из левого верхнего угла расчетной области, которая взаимодействует впоследствии с отраженной волной. Аналогичная рис. 10 картина изолиний плотности для задачи 2 получена в работе [17] для монотонной схемы с экстраполяцией вдоль характеристик. Можно предположить, что формирование вторичной нефизической волны связано с несогласованностью конечно-разностной схемы с сильным разрывом в краевых условиях. Однако такие же возмущения замечены при расчете сверхзвукового обтекания цилиндра, где разрыв образуется в центральной части расчетной области, а значит, не могут быть объяснены влиянием граничных условий. Таким образом, нефизический "провал" является чисто схемным эффектом, характерным для всех рассмотренных нами схем.

Истории сходимости численных алгоритмов при решении задачи 2 представлены на рис. 11.

Рис. 11. Сходимость TVD-схем на задаче 2

= 30°).

Рис. 12. Сходимость TVD-схем на задаче 2 (^ = 29°).

Как и следовало ожидать, усложнение задачи по сравнению с задачей 1 приводит к ухудшению сходимости всех схем. При этом сходимость методов HarMIN и D2B1, имевших наилучшие результаты на задаче 1, замедляется примерно в 4 раза. Заметим, что симметричная схема SYM, сравнимая с ними по точности, не позволяет достичь сходимости даже Err = 10-4. Невязка схем высокого класса точности, таких как D2B3, а также алгоритма HarVL не падает ниже значения Err ~ 10-6, чаще всего выходя на периодический режим колебаний.

Для иллюстрации чувствительности схем к малым изменениям входных параметров, таких как расположение точной линии разрыва по отношению к узлам сетки, был проведен расчет задачи 2 при угле наклона скачка ^ = 29° (это достигается путем соответствующей модификации граничных условий на верхней входной границе). Результаты для этого случая представлены на рис. 12. Сравнение кривых, соответствующих схемам HarSB, UPYee, OHarW1 на рис. 11 и рис. 12, обнаруживает существенную зависимость их сходимости от небольших изменений входных данных, причем без какой либо закономерности. Среди представленных на рис. 11 и 12 алгоритмов выделяется схема Чакравати с ограничителем Ван Лира (D2VL), которая в обоих случаях позволяет достичь сходимости

невязки до 10-14 менее чем за 400 итераций. Таким образом, имея вполне приемлемую точность, D2VL превосходит все остальные рассмотренные TVD-схемы в плане сходимости. Это явление нельзя объяснить одной лишь гладкостью ограничителя (18) в точке 0j = 1 (см. рис. 1), так как схема HarVL также использует этот ограничитель.

Заключение

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

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

2. Высокоточные TVD-схемы, такие как схема Хартена с искусственным сжатием при и = 2, и схема Чакравати третьего порядка при в = 3, вообще говоря, не дают установившегося решения. В тестовых задачах предельная невязка решения стационарных разностных уравнений Err ~ 10-4.

3. Среди рассмотренных схем лучше всего сочетает в себе достаточную точность воспроизведения разрывов и высокую скорость сходимости к стационарному решению TVD-схема Чакравати с ограничителями Ван Лира.

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

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

[1] Стрелец М. Х., Шур М. Л. Метод масштабирования сжимаемости для расчета стационарных течений вязкого газа при произвольных числах Маха. ЖВМ и МФ, 28, 1988, 254-266.

[2] Choi D., MERKLE C. L. Application of Time-Iterative Scheme to Incompressible Flows. AIAA J., 23, 1985, 1518-1524.

[3] Harten A. High Resolution Schemes for Hyperbolic Conservation Laws. J. Comput. Phys, 49, 1983, 357-393.

[4] Yee H.C., Warming R. F., Harten A. Implicit Total Variation Diminishing (TVD) Schemes for Steady-State Calculations. Ibid., 57, 1985, 327-360.

[5] Yee H. C. Construction of Explicit and Implicit Symmetric TVD Schemes and Their Applications. Ibid., 68, 1987, 151-179.

[6] Moon Y. J., Yee H.C. Numerical Simulation by TVD Schemes of Complex Shock Reflections From Airfoils at High Angle of Attack. AIAA Paper 87-0350, 1987.

[7] Белецкий Ю. М., Войнович П. А., Ильин С. А. и др. Сравнение некоторых квазимонотонных разностных схем сквозного счета. 1. Стационарные течения. Препр. ФТИ им. А. Ф. Иоффе, №1383, Л., 1989.

[8] Ильин С. А., Тимофеев е. В. Сравнение некоторых квазимонотонных разностных схем сквозного счета. 2. Линейный перенос возмущений. Там же, №1550, Л., 1991.

[9] Sundar K., Sriramulu v., Ramakrishna M. Comparison of High Accuracy TVD Schemes to Flows Containing Strong Shocks. AIAA J., 33, № 11, 1995, 2087-2091.

[10] Takakura Y. ET AL. On TVD Schemes for 3D Euler Equations in General Coordinates. Int. J. for Numerical Methods in Fluids, No. 9, 1989, 1011-1024.

[11] Roe P. L. Approximate Riemann Solvers, Parameter Vectors and Difference Schemes. J. Comput. Phys., 43, 1981, 337-372.

[12] Chakravarthy S.R., Osher S. A New Class of High Resolution TVD Schemes for Hyperbolic Conservation Laws. AIAA Paper 85-0363, 1985.

[13] Чакравати С. Р., Жем К.-Й. Расчет трехмерных сверхзвуковых течений с дозвуковыми зонами на основе уравнений Эйлера. Аэрокосмическая техника, №11, 1987, 22-35.

[14] Shi J., TORO E. F. Fully Discrete High-Order Shock-Capturing Numerical Schemes. Int. J. for Numerical Methods in Fluids, 23, 1996, 241-269.

[15] Van Leer B. Towards the Ultimate Conservative Scheme. 2. Monotonicity and Conservation in a Second-Order Scheme. J. Comput. Phys., 14, No. 2, 1974, 361-376.

[16] Грязин Ю.А., Черный С. Г., Шаров С. В., Шашкин П. А. Об одном методе численного расчета решения трехмерных задач динамики несжимаемой жидкости. Докл. РАН, 353, №4, 1997, 478-483.

[17] Lacor C., Hirsch Ch. Genuinely Upwind Algorithms for the Multidimensional Euler Equations. AIAA J., 30, No. 1, 1992, 56-63.

Поступила в редакцию 7 января 2000 г.