CC BY
47
Е.И. Васильев, A.C. Демин, 2006
ПРИКЛАДНАЯ МАТЕМАТИКА
УДК 519.6 : 533.7
3D ВАРИАНТ W-МОДИФИКАЦИИ МЕТОДА ГОДУНОВА ДЛЯ НЕСТАЦИОНАРНЫХ ТЕЧЕНИЙ СОВЕРШЕННОГО ГАЗА
Е,И. Васильев, A.C. Демин
Предлагается модификация метода Годунова повышенной точности для трехмерных нестационарных течений совершенного газа. Доказано свойство позитивности предложенной схемы для линейного скалярного уравнения переноса. Проведено тестирование метода на примере задачи о запуске трехмерного профилированного сопла квадратного сечения.
Введение
В последнее время в вычислительной газодинамике широкое распространение получили монотонные явные методы повышенной точности, основанные на существенно нелинейных разностных схемах. Основная идея подобных схем заключается в консервативной поправке потоков через границы ячеек для линейной схемы первого порядка [1, 4, 5, 6]. Нелинейный характер поправок, а также переменный шаблон их вычисления позволяют обойти известную теорему Годунова [2] и совместить второй порядок точности схемы с монотонностью.
В данной работе предлагается модификация метода Годунова с вычислением поправок на трехмерном переменном шаблоне.
Пусть гиперболическая система уравнений записана в дивергентном виде:
А
угь + РЫя + + Н(\у)2 = О,
где \у(х, у, г, ¿) — вектор из т компонент, Р(\у), 0(\у), Н(\у) - вектор-функции потоков. Тогда метод Годунова 1-го порядка [3], примененный к модифицированному уравнению
+ Р(\У + а)ж + + /3)у + Н(лу + 7)г = О,
даст приближенное решение исходного уравнения со вторым порядком точности по пространству и времени, если
Ах . At
<* = -^-sign(Fw)wI + -у w*,
_ Ay . At
P = -ysign(Gw)wJ/ +
Az . At
7 = —sign(Hw)wz 4- — wt,
где Fw, Gw и Hw — матрицы производных.
В последних формулах производные по времени wt можно заменить на пространственные производные с помощью исходного уравнения:
а = ^sign(Fw)wI - ^F(w)B - ^G(w)у - ^H(w)2,
/3 = ^sign(Gw)wy - yF(w)x - yG(w)y - ^H(w)2, * (1)
7 = ^sign(Hw)wz - ^F(w)e - ^Gfw)^ -
Аппроксимация величин a, /3 и 7 при численном решении осуществлялась нелинейным образом на трехмерном W-шаблоне, ориентация которого зависит от направления скорости потока (собственных чисел матриц Fw, Gw и Hw).
1. Анализ схемы для линейного скалярного уравнения
Рассмотрим трехмерное скалярное уравнение переноса с постоянными коэффициентами:
dw .dw dw dw . А ...
— + A— + ц— + = 0, A, ¡x, rj = const. (2)
Пусть сетка выбрана равномерной и прямоугольной, причем направление роста индексов ячеек (j,i,k) совпадает с направлением осей координат (x,y,z). Схема Годунова в случае положительных коэффициентов уравнения (2) имеет вид:
IVjik Wjik , \ i Wjik Wjti— , „Wjik 1Uji,k—1 n
------------г Л----------------г U-----7--------Г Tl------7 — U,
At Ax Ay Az
где Wjik и Wjik — значения искомой функции w(x,y,z) в ячейке (j,i,k) в момент времени t и t + At.
Для более общего случая введем обозначение для индексов, «сдвинутых назад» — j_, г_, /г_, и индексов, «сдвинутых вперед» по потоку — j+, г+, к+:
j- =j~ sign (A), i- = i-sign(/i), к— — к sign (77),
j+ = j + sign(A), г+ = г + sign(^), k+ = к + sign(ij).
Для каждой ячейки с номером (j, г, к) рассмотрим разности против потока вдоль осей х,у a z - pjik, qjik, rjik:
Pjik Wjik 'UJj-iki
Qjik Wjik Wji__ki (4)
Tjik ~ Wjik Wjifc_.
В этих обозначениях разностная схема Годунова (3) для (2) запишется:
Wjik — Wjik VxPjik VyQjik Vz^jiki
(5)
где
Схема (5) устойчива и позитивна при выполнении следующих условий:
Vx > о, Vy > 0, vz > О, V = vx + vy + vz < 1
(6)
где V — сумма чисел Куранта для различных направлений.
Введем величины, аппроксимирующие выражения (1) с помощью разности против потока:
Модифицированную схему запишем в виде, аналогичном схеме Годунова (5):
где разности pjik, qjik и вычисляются аналогично (4), но с поправками:
Поправки dxwjik, dyWjik и dzWjik определяются как нелинейные средние из пар величин (7) для ячейки (j,i,k) и ячейки, находящейся впереди по потоку:
(F'lVjik = mid (Oijiki Qj+ik)) dyWjik — mid (/3j, Pji+k)i 9 Wjik mid [pfjiki (Ю)
mid (a, b) — некоторая нелинейная функция осреднения а и b.
Отметим, что величины otjik,aj+ik, /3jik,Pji+k и 7jifc>7jifc+ вычисляются на шаблоне, зависящем от знаков А, ц и rj, с ориентацией ветвей против потока. На рис. 1 изображены два варианта упомянутого шаблона при различных положениях вектора скорости переноса (A, /lx, 77).
Обоснование позитивности схемы дает следующая теорема.
Теорема 1. Если функция mid (a, b) такова, что для любых а, Ъ, с
®jik = (р~ vхР ~ vyq ~ vzr)jik,
Pjik = (q-~ vyQ - vxP ~ yzr)jik,
Ijik = (r - Vzr - uxp - uyq)jik.
(7)
Wjik — Wjik VxPjik VyQjik VzTjik,
(8)
Pjik = (w + dxw)jik ~{w + dxw)j_ik, Qjik = (w + dyw)jik - (w + dyw)ji_k, rjik = (w + dzw)jik ~ {w + dzw)jik_.
(9)
|mid(a,b) — mid (c,a)| < r\a\,
(11)
2
(12)
r =
l + v
то при условии (6) разностная схема (7)-(10) позитивна.
Рис. 1. Переменный шаблон для вычисления поправок
Доказательство. После подстановки (9) в схему (8) имеем:
- 'ujjik = Wjik vx(,Pjik d Wjik д Wj_ik)
L'yiQjik ~t“ dyWji_ k ) Vz(j"jik + d Wjik d 'UJjik- ) •
Для разности поправок из (10) и (11), с учетом того, что otjik = а^_)+ік,
ßjik ~ ßj(i~)+k И Tjifc = ')ji(k-)+1 ПОЛуЧИМ.
IdxWjik - dxWj_ikI < r\ajik\, Idywjik - dyWji_k\ < r\ßjik\, \dzwjik - dzwjik_ \ < r\^jik\. Последнее можно переписать в виде: dxwjik - dxwj_ik = irajik) dywjik - dywji_k = (rßjlk, dzwjik - dzwjik_ = ^njik, где С, С и ф — некоторые числа:
-1<£<1, -1 <С< 1, (13)
Тогда схему (8) можно представить в виде, аналогичном схеме Годунова (5):
lUjik — jik VxPjik VyQjik ^z'fjik)
где
% = vx[ 1 + Ur( 1 - ux) - kruy - Ыгиг],
Vy = Uy[l + kr{l - Vy) - krux - Uruz],
Vz = ug[l + \фг( 1 - uz) - \t,rvx - \CrUy].
Из последних формул с учетом (13) получаем следующие оценки:
"х > Vx[l - Ml -Vx + Vy + I/,)],
Vy>Vy\l-\r(l-Vy + Vx + Vz)\, П4,
Vz > vz[l - |r(l - Vz + Vx + Vy)],
l-Vx-Vy-Ug>(l-Vx-l'y- UZ)[1 - l{vx + Vy + Uz)],
из которых, используя (6) и (12), следуют условия, аналогичные условиям (6):
их > 2их[^^\ >0, Уу > 2г/у[-^;] >0, г/2 > 2^г[^^;] > 0,
1 - дх - ду - > (1 - г/)^] > 0.
Следовательно, схема (7)—(10) является позитивной, как и исходная схема Годунова.*
Выбор функции г вида (12) для доказательства теоремы не является оптимальным. Рассмотрим выражения (14), полученные в ходе доказательства теоремы, и
потребуем, чтобы их > 0, иу > 0, иг > 0 и 1 — их — иу — Рг > 0. Тогда, выражая
из неравенств г, получим три не поглощающих друг друга неравенства:
г <тх = ----------------, r<rv~---------------------, r<rz = ------------------—.
1 1SX “I" ISy Vz 1 l/y -f- l/x -f- Uz 1 1SZ “Ь Vx “f" Uy
(15)
Таким образом,-в качестве функции г можно взять
г - min(rx,ry,rz).
Это обеспечит выполнение условий (15) и расширит диапазон значений г. Обеспечить более мягкие ограничения на г можно лишь за счет рассмотрения векторной функции г = (тХ1 Гу, rz), где компоненты гх, гу и rz в общем случае различны.
Приведем несколько функций [1], удовлетворяющих условиям (11), (12):
mid (а, Ь) = mms(a,b) = ( min(\a\, |6|)sign(a), если ab > 0, v ' 4 ' [0, если ab < 0.
.,/ , / ,\ r(\ab\ + ab)sign(a)
mid (a, b) = mhsr(a, 6) =-------— ■ = = (17)
|а + Ь\ + уДа + Ь)2 — 4аЬг(2 — г)
При г = 1 последняя функция обращается в mms(as Ь). При г = 2 она обращается в функцию «гармонического среднего», которую для одномерных задач впервые предложено использовать в [6].
Условия (11), (12) для последних двух функций (16), (17) выполняются, так как справедливо соотношение:
\mhsr(a,b)\ < r\mms(a,b)\.
А
2. Вычисление поправок для векторного уравнения
Рассмотрим гиперболическую систему уравнений с тремя пространственными переменными, записанную в недивергентной форме:
dw ^dxv dw
Ж + х & + у Ц + zTz = °’
где w — вектор из т компонент, а X, Y и Z — матрицы, зависящие от w. Собственные числа упомянутых матриц действительны, а множество собственных векторов образует базис, то есть имеет место разложение:
X = RxLRxl, Y = RYMRy\ Z = RzNRTzx,
где Rx, Ry, Rz — матрицы собственных векторов, a L, М и N — диагональные матрицы собственных чисел (А/), (щ) и (щ), I — 1,2
Пусть сетка выбрана равномерной и прямоугольной, причем направление роста индексов ячеек (j,i,k) совпадает с направлением осей координат (x,y,z). Аналогично скалярному случаю введем обозначение для индексов, «сдвинутых назад» - ji, ц, ки и индексов, «сдвинутых вперед» по потоку — jp, ip, кр:
ji = 3 ~ sign(Ai), ii = i - slgn(m), kt = k-sign(rji), jp = j + sign(Ap), ip = i + s\gn(vP), kp = к + signup).
Рассмотрим для каждой ячейки с номером (j,i,k) векторы характеристических разностей против потока вдоль осей х, у и z — pjik, qjik и ([ ]/ — обозначение l-й компоненты вектора):
[рjik\i = [RxHvfjik - wjlik)]i, hjikh = [Ry^Wjik - wj4|fc)]i,
[rjik\i = [R^iyvjik ~ W
где — вектор параметров в ячейке (j,i,k) в момент t.
Введем величины, которые с точностью до матричного сомножителя аппроксимируют выражения (1) с использованием разностей против потока:
К*]i = [р¡¡к - ШЦрщ, - ЩЯх'Яг\М]Чцк - ^R-xlRz\N\rm}h
[/3jtt]l = - ШМ\Ч1Л - £R?Rx\L\pjik - &R?Rz\N\riik]„ (18)
Ь,Л = [fjtt - Sl^|r№ - - gHi1 Яу|М|%й]„
И окончательно для каждой ячейки находим три вектора поправок dxwjik, dywjik, dzv/jik, как нелинейные средние из пар величин (18) для ячейки (j,i,k) и ячейки, находящейся впереди по потоку:
dxwjik = \Rx- mid {a.jik, Oijpik),
dywjik = |Ry ■ mid ((3jik, f3jipk), (19)
= f • mid hjiknjikp)-
Нелинейная функция осреднения mid (a, b) для векторных величин применяется покомпонентно. Значения коэффициентов матриц вычисляются по параметрам той ячейки, в которой производится вычисление поправок. Поправки постоянны для данной ячейки и добавляются к основным параметрам перед вычислением потоков через соответствующие границы ячейки. В методе Годунова для вычисления потоков используется решение задачи Римана.
3. Схема численного решения трехмерных уравнений газодинамики
Дивергентная форма системы дифференциальных уравнений, описывающих трех мерное нестационарное течение невязкого совершенного газа, имеет вид:
<эи эт(и) №(и) ан(и)
ді
+
дх
+
ду
+
дг
= 0,
и =
/ р ^ ( Ри ^ ( рп ^ / рги \
ри ри2 + р рии рии>
ру , р = риу , о = рь2 + р , Н = рьи>
ри] риио рут рго2 + р
V е ) \ (е + р)и ) \ (е + р)^ ) \ (е + р)и) у
(20)
(21)
где р и р — плотность и давление соответственно, и, V, ги — компоненты вектора скорости в трехмерном пространстве, е = р/(7 — 1) + 0,5р(и2 + V2 + ги2) — полная энергия единицы объема, 7 — показатель адиабаты. Вектор-функции Р(и), 0(11) и Н(и) называются функциями потоков.
Рассмотрим отдельную подвижную ячейку 7} в течение времени от ¿п до <п+1 — %г + ^ как объемное тело V,- в четырехмерном пространстве ОхугЬ. Нижнее и верхнее основания четырехмерного тела Г" и Г"+1 представляют собой ячейку 7} в различные моменты времени. Считаем, что при своем движении ячейка остается шестигранником, причем вершины в течение шага Д£ двигаются с постоянной скоростью (рис. 2).
* Рис. 2. Подвижная ячейка в пространстве (ж, у, г, <)
Интегральный аналог уравнений (20), примененный к этому телу, можно записать в следующем виде:
® (\Jrrit + Ттх + Ъту + Нт2)дУ = 0. (22)
Щ
Слева интеграл по поверхности четырехмерного тела = Т? и Г"+1 У(У Т^к),
к
где Т^, /с = 0,5 — боковые трехмерные грани, (mt,mx,my,mz) — компоненты внешней нормали ш к поверхности этого тела. С учетом того, что Т” и Т"+1 параллельны гиперплоскости Охуг, формулу (22) можно записать в виде:
/иж + j иmtds + ^ (игл* + ¥тх + вту + Шт2)дУ = 0.
Ч т;+1 иък
к
Заменяя на каждой из граней подынтегральные величины на постоянные (например, считая их равными значениям в центрах граней, то есть в точках пересечения средних сечений), можно записать следующее приближенное равенство:
и'п+1Уп+1 = и пуп _ , * (23)
к
где .
= И]к ! гщбУ + ¥:1к ! тхдУ + О,, 'I гпуёУ + П]к J тгвУ. (24)
Т3к Тп: Т3к Т3к
В (23) V-1 и У™+1 — объемы верхнего и нижнего оснований. Величины 0^ представляют собой векторы потоков через jk грань ячейки за промежуток времени АЬ. Интегралы в (24) представляют собой алгебраические значения (то есть с учетом знака) объемов проекций к-го шестигранника на соответствующие подпространства Охуг, Оху1, Охг1 и ОухЬ. Далее индекс ¿к временно опускаем для краткости.
На рис. 3 изображена подвижная грань ячейки в моменты t к t + АЬ (затененные четырехугольники), вер — среднее сечение шестигранника, 5 — площадь
среднего сечения, п = (пх,пу,пг) — единичный вектор нормали к среднему сечению, направленный вдоль проекции вектора гп на гиперплоскость t = £га + Д£/2, —
объем проекции шестигранника на гиперплоскость Охуг, d — осредненная скорость движения грани ячейки.
Через векторы а, Ь и с, которые соединяют середины противоположных сторон подвижной грани ячейки, можно выразить все ранее упомянутые параметры:
Б = \а хЦ, (пх,пу,пг) = -(о х Ъ), П =<а,Ь,с>, (25)
а также выразить интегралы:
J mtds = Г2 = —SAtd, j тхдУ = ЗпхАЬ,
Т3к
J mydV = SriyAt, J mzdV = SnzAt.
Tik T3k
Рис. 3. Подвижная грань ячейки
После подстановки в (24) получим:
О
■Зк
6^. А.і(Тпх + (лпу + Н пг — и в)
Зк'
С учетом (21) последнее выражение можно записать более подробно:
Зк
( р{ип -6) \
ри(ип - (1) + пхр
ру(ип - й) + ПуР
ри)(ип - (Г) + щр \ ре{ип - д) + ипр
ип — п-хЫ + пуУ + пгт.
(26)
Зк
Для вычисления потоков по (25)—(26) необходимо указать способ определения величин (р,р,и,ь,ги)зк на каждой к-й грани ¿-й ячейки. В методах типа Годунова эти величины определяются как решение вспомогательной задачи о распаде разрыва (задача Римана) на грани с нормалью п и движущейся со скоростью с1. В методе Годунова 1-го порядка в качестве исходных данных для задачи Римана используются параметры в соседних ячейках с временного слоя £ = £п. В '^модификации метода Годунова к этим параметрам перед решением задачи Римана добавляются корректирующие “поправки, которые обеспечивают повышение порядка метода до 2-го по пространству и времени.
Рассмотрим более подробно алгоритм вычисления поправок для уравнений движения совершенного газа. С целью экономии вычислений используем недивергентную форму уравнений Эйлера, которые в декартовых координатах имеют вид:
Ж + ХЖ + УЖ+ Ж = 0'
(27)
w =
Y =
V
о
V
О
о~1
О
и
V
W
О
V w
( и О О О О
О 7р
V
О
О
о
V
О
о \ о о о v}
и
ог1
0
О
р
IV
и
О
О
О
О
о \ о
о о
и О
о
z =
/ w О О
о
Vo
о
w
О
о
„-1
и j
О о
о
w
О
О
О
о
w
О
р \
1Р
о о
w j
(28)
Система уравнений (27) используется при вычислении поправок по формулам (18)—(19).
Процедура вычисления поправок, описанная в главе 2, не использует ортогональности сетки или координат (а только равномерность) и не изменится, если х, у и z — криволинейные координаты, а сетка, построенная из координатных линий, является равномерной именно в этих криволинейных координатах.
В изначальном уравнении (27) при переходе в криволинейные координаты (£> Vi О- которые связаны с сеткой, произойдет преобразование матриц системы уравнений (28) и соответствующее этому изменение в процедуре вычисления поправок.
Пусть у нас есть некоторая криволинейная неравномерная и неортогональная подвижная сетка с четырехугольными ячейками. Обозначим через (£,??, С) криволинейные координаты, в которых сетка неподвижна, равномерна, и направление роста £, r¡ и ( совпадает с направлением роста индексов ячеек. Координатными линиями этой криволинейной системы координат являются сеточные линии, образующие сетку. Предположим наличие дифференцируемой зависимости:
х = ж(£,?7, С.*). z =
(29)
Введем единичные векторы а = (а1} а2, аз)т> Р = (Pi, @2, Рз)т и 7 = (7ь 72,7з)Г. ортогональные координатным поверхностям £ = const, т? = const и £ = const соответственно, и вектор скорости d перемещения криволинейной системы по отношению к неподвижным декартовым координатам:
УчЧ - ЧУ< ЧУс - УаЧ Уе2п ЧУ у
а ~ а а Zr/X£ — XVZ£ - ß = ~ß чч - чч , 7=7 z^xr¡ — x$zv
ЗД Уг]Х^ щх^ - х^_ X&r, - VtX
а
d =
xt
Vt
Zt
\/Т^£С - ^У<)2 + (^а?с ~ хуЧ)2 + (хуУ< ~
Р = \Лчус ~ у^ч)2 + (чч - чч)2 + (у&с ~ ЧУ<д21
7 = у/Лу&, - г^Уг,)2 + (г^Хг, - х^)2 + {х^уп - у^)2,
С учетом приведенных выше обозначений уравнение (27) преобразуем к следующему виду:
дуг . 9w ^Эуг „дч?
т+Ащ+в^ + сЖ = 0'
г
где
А — — (оііХ + о:2У + — {{ос • <£)Е)),
и
в = ^(дх + &у + дг - ((3 ■ 3)Е)),
С = 1(ЪХ + ъУ + 7з2 - ((7 • 3)£0),
2/77
и якобиан преобразования У =
ХС 2С
В разложении матриц Л, В, С,
Л = ВдАаВд1, В = ВвАвВ^, С = ВсЛсВ^, матрицы собственных векторов и собственных чисел устроены схожим образом:
1 Пі тг2 1
с2 0 0 0 с2
Ядг = -Пі с/р 0 —тг3 П-2 Пі с/р
-п2с/р о со е -пх П2с/р
-щс/р -гг2 пі 0 и3с/р
'о 1/2 с2 - -пір/2с -н2р/2 с —пзр/'2с
П\ -Пі/с2 0 п3 — «2
II п2 -п2/с2 -Пз 0 нг 5
Щ -713/с2 п2 -Пі 0
0 1/2 с2 пір/2с п2р/2с п3р/2с _
(д-гі) — с 0 0 0 0
0 й • (д — с?) 0 0 0
0 0 и - (д- 3) 0 0
0 0 0 п ■ (д-сО 0
0 0 0 0 и - (д — (І) + с_
Ад' =
Здесь индекс ІУ принимает одно значение из (Л, Б, С), а вектор и = (пь п2, п3) соответственно принимает значение а, ¡3 или 7. Вектор д = (и,ь,и>)т — вектор скорости газа, с - скорость звука в газе (с = утр/р).
Структура произведений матриц В^Вв, В^Вд, В^Вс, В^1 Ва, Вд1 Вс и входящих в формулы для поправок, довольно проста. Так, например, имеем: *
Вд1Вв = Р-(Э, В~в1ВА = Р + (Э,
(1 + в)/2 ООО (1-з)/2 О 5 О О О
О 0 я О О
О 0 0 5 0
(1-а)/2 ООО (1 + в)/2
0 Р^2з/2с -pyiz/2c pzW2c 0
~су23/р 0 —у 12 —^13 СУю/р
CVu/p У\2 0 — ^23 -CV13/P
-СУ12/Р ¿73 ^23 0 СУ12./р
0 ~ Р^2з/2с ¿и/13/2с -pz/12/2c 0
где
5 = ОііРі + 0-2^2 + «зДз) /зд\
^12 = «і/?2 — Оі2/Зі, ¡73 — а\Рг — СИз/?і, У22, — &2Р3 ~ »з/^-
Для получения произведений К^Яс и Д^:Яв достаточно в выражениях (ЗО) заменить а на (3, а /3 на 7.
Для практических вычислений удобно принять, что в криволинейных координатах (£,77,С) узлы сетки имеют целочисленные координаты. Это правомерно, так как преобразование (29) не единственно и определяется с точностью до масштаба измерений вдоль (£,??, С)- Тогда, переходя к дискретному аналогу ячейки как шестиграннику (см. рис. 4), получим, что а, /3 и 7 - единичные векторы, перпендикулярные к средним сечениям ячейки, а, (3 и 7 — площади средних сечений ячейки, J — объем ячейки. '
Рис. 4. Трехмерная расчетная ячейка
4. Результаты тестирования метода
Работоспособность метода рассмотрим на примере задачи о формировании течения в профилированном трехмерном канале в режиме трубы Людвига, то есть, когда диафрагма, обеспечивающая начальный перепад давления, располагается около выходного сеч‘ения канала. Геометрическая конфигурация канала изображена на рис, 5. Рассматриваемая область канала имеет длину Ь = 2,5, диафрагма, обозначенная к отмеченная на рисунке жирной линией, располагается на расстоянии 0,1Ь от выходного сечения и обеспечивает начальный перепад давления, равный 20. Температура по разные стороны от диафрагмы одинаковая.
В рассматриваемой геометрии сопла участок расширения заменен параллельным каналом постоянного сечения. Это сделано для того, чтобы устранить образо-
вание вторичных скачков в этом участке. Реальный порядок точности метода определяется лишь на гладких течениях, а вторичные скачки не позволяют это сделать достаточно аккуратно.
А
Рис. 5. Конфигурация сопла в тестовой задаче
В начальный момент газ по обе стороны диафрагмы покоится. После мгновенного разрыва диафрагмы образуется ударная волна и контактный разрыв, которые распространяются вправо и достаточно быстро выносятся из канала. Влево внутрь трехмерного канала распространяется нестационарная волна разрежения, за счет которой осуществляется запуск сопла. В расчетах на стенках канала использовалось условие непротекания.
На рис. 6 изображены изолинии плотности в плоскости симметрии канала на момент времени £ = 2 после разрыва диафрагмы. К этому моменту времени волна разрежения прошла переходной участок от широкого канала к узкому. Около стенки на входе в узкий параллельный канал наблюдаются повышенные градиенты плотности. Влияние трехмерных эффектов в этой области можно проследить по рис. 7, где на тот же момент времени изображены изолинии плотности в поперечных сечениях канала с координатами х = 1,-5, х = 1,6 и 1 = 1,7 с одинаковым шагом Ар. Сечение а: — 1, 6 соответствует точке сопряжения, то есть началу параллельного участка. Хорошо видно, что в этом сечении градиенты в угловой области на стыке двух стенок в несколько раз больше, чем градиенты в плоскости симметрии, что говорит о ярко выраженном трехмерном характере течения. Вправо по потоку неоднородность течения быстро уменьшается, в сечении х = 1,7 высокие градиенты отсутствуют, ь На данный момент времени течение внутри канала непрерывно и энтропия теоретически должна быть равна начальному значению, то есть величина 5 = р/р1 должна быть равна начальному 50 = Ро/Ро ■ Этот факт можно использовать для анализа точности метода путем сравнения численного и теоретического значения энтропии.
Рис. 6. Изолинии плотности в контрольный момент времени
Рис. 7. Изолинии плотности с шагом Ар — 0,12 в поперечных сечениях сопла х = 1,5, х = 1,6 и х = 1,7 в контрольный момент времени
В расчетах использовалась неподвижная сетка, включающая N ячеек по длине сопла и N/2 по ширине и высоте. На рис. 8 для различных значений N = 30 4-180 изображено распределение ошибки AS/h2 вдоль средней линии в момент t = 2.
Ошибка в энтропии AS = (S — S0/Sq) отнесена к квадрату горизонтального шага сетки h = L/N. Хорошо видно, что, начиная с N — 90, поведение кривых меняется незначительно. Это означает, что ошибка в энтропии AS ~ Ch2, где С не зависит от N.
Результаты рис. 8 объединены в один график на рис. 9а, где представлена зависимость нормированной средней ошибки в энтропии от числа ячеек.
Видно, что при N > 90 рост ошибки замедляется. На рис. 96 изображена интегральная ошибка в энтропии, нормированная на h1'8. Из сравнения кривых на рис. 9 можно сделать вывод, что порядок точности метода лежит в промежутке 1,8-2.
При изменении h пропорционально изменяется и шаг по времени Ai. Следовательно, учитывая существенную нестационарность течения, близкий ко второму порядок точности имеет место не только по пространству, но и по времени.
Рис. 8. Распределение нормированной ошибки в энтропии вдоль средней линии сопла на различных сетках
Рис. 9. Зависимость интегральной нормированной ошибки энтропии на средней линии сопла от количества ячеек
В заключение отметим вычислительную эффективность предложенной схемы, которая связана с двумя основными моментами.
1. Отсутствие расщепления по пространственным координатам значительно снижает суммарный объем пересылаемых данных при параллельных расчетах задач на нескольких компьютерах.
2. Использование неконсервативной формы уравнений для вычисления поправок приводит к матрицам простой структуры с большим количеством нулевых элементов, что позволяет сделать довольно экономичный алгоритм вычисления поправок.
Summary
3D VARIANT OF THE W-MODIFICATION OF THE GODUNOV METHOD IN NON-STATIONARY FLOWS OF PERFECT GAS
E.I. Vasilev, A.S. Demin
A modification of the Godunov method with increased accuracy in three-dimensional flows of a perfect gas was suggested. It has been proved the positivity property of the suggested scheme for a linear scalar transport equation. Testing of the method has been made on the problem of the launching of three-dimensional grading nozzle with square cross-section.
Список литературы
1. Васильев Е.И. '\У-модификация метода С.К. Годунова и ее применение для двумерных нестационарных течений запыленного газа // Журнал вычислительной механики и математической физики. 1996. Т. 36. № 1. С. 122-135.
2. Годунов С.К. Разностный метод численного расчета разрывных решений уравнений гидродинамики // Матем. сб. 1959. Т. 47. С. 271-306.
3. Численное решение многомерных задач газовой динамики / Годунов С.К., За-
бродин А.В., Иванов М.Я. и др. М.: Наука, 1976. 400 с.
4. Colella P. Multidimensional Upwind Methods for Hyperbolic Conservation Laws
// J. Comput. Phys. 1990. V. 87. P. 171-200.
5. Harten A. High Resolution Schemes for Hyperbolic Conservation Laws
// J. Comput. Phys. 1983. V. 49. P. 357-393.
6. Leer B. van Towards the ultimate conservative difference scheme. V. A second-order
sequel to Godunov's method // J. Comput. phys. 1979. V. 32. P. 101-136.
