CC BY
148
Вычислительные технологии
Том 12, № 3, 2007
МОДИФИКАЦИИ АЛГОРИТМОВ РАСЩЕПЛЕНИЯ ДЛЯ РЕШЕНИЯ УРАВНЕНИЙ ГАЗОВОЙ ДИНАМИКИ И НАВЬЕ-СТОКСА*
В.М. Ковеня, А. Ю. Слюняев Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: [email protected], [email protected]
Various finite difference schemes with the effective splitting of operators aimed at the numerical solution of the gas dynamics and the Navier—Stokes equations applied for a viscous compressed heat conducting gas flow are offered in this paper. The generalization of schemes for the case of the curvilinear coordinates is given. Results of the numerical modeling of the internal flows in a channel with the shock waves formation as well as the shock-to-shock and shock-to-wall interaction and the formation of separation zones along with other features of the flow are presented.
Введение
Уравнения газовой динамики и Навье—Стокса вязкого сжимаемого теплопроводного газа служат базовой моделью для широкого класса задач аэро- и гидродинамики. Их решения содержат, как правило, узкие зоны больших градиентов и другие особенности, такие как пограничные слои, висячие скачки, отрывные зоны и т. д. Поэтому их численное решение представляет значительные трудности, что вызывает необходимость разработки и использования специальных численных алгоритмов. Эти алгоритмы должны обладать достаточной точностью, большим запасом устойчивости, удовлетворять свойствам консервативности и экономичности с возможностью получения решения задачи за разумное время на существующих ЭВМ. Использование явных разностных схем для численного решения уравнений Навье—Стокса и, в ряде случаев, уравнений газовой динамики неэффективно в силу жестких ограничений на соотношение временного и пространственных шагов расчетной сетки, особенно при нахождении стационарного решения методом установления, поэтому часто используются неявные разностные схемы, свободные от ограничений на устойчивость или имеющие более слабые ограничения. Обзор наиболее употребительных разностных схем приведен, например, в работах [1-7].
При построении и реализации неявных разностных схем наиболее часто используются методы факторизации и расщепления (см. [1-8]), которые позволяют свести решение исходных многомерных задач к последовательному (или параллельному) решению их одномерных аналогов. Однако даже решение одномерных задач по неявным разностным
* Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (грант № 05-01-00146).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.
схемам приводит к необходимости обращения матриц, размерность которых возрастает с ростом числа уравнений или с увеличением размерности задач по пространству. Поэтому остается актуальным построение численных алгоритмов — их реализация экономична и сводится, к примеру, к скалярным прогонкам или к "схеме бегущего счета" для каждого пространственного направления, как, например, в схемах расщепления [3, 8]. Известно (см. [8, 10]), что введение расщепления или факторизация операторов в исходной многомерной задаче приводит к появлению дополнительных членов в разностной схеме — диссипатив-ных членов и членов более высокого порядка, и, как следствие, к ухудшению свойств численного алгоритма. Это и является платой за экономичность метода. Поэтому расщепление операторов следует выбирать таким образом, чтобы минимизировать влияние этих членов.
В настоящей работе предлагаются и исследуются разностные схемы, основанные на расщеплении исходной многомерной задачи на их одномерные аналоги с последующим расщеплением одномерных задач таким образом, чтобы реализация схем на дробных шагах сводилась к скалярным прогонкам, т. е. они были бы экономичны по числу операций на узел сетки, построенные схемы были бы безусловно устойчивыми или имели бы слабые ограничения на устойчивость, а влияние расщепления было бы минимальным, т. е. их свойства были бы близки к свойствам нефакторизованных схем. Некоторые схемы минимальной диссипации для решения одномерных задач газовой динамики рассматривались ранее в работе [10]. Ниже предлагаются и исследуются разностные схемы для многомерных уравнений газовой динамики и уравнений Навье—Стокса сжимаемого теплопроводного газа в декартовых и криволинейных координатах. Во второй части работы приведены результаты тестирования предложенных алгоритмов и результаты расчетов течений вязкого газа в плоском канале, в том числе с учетом вдува газа с части его поверхности.
1. Разностные схемы для уравнений газовой динамики
Построение экономичных разностных схем, основанных на идеологии расщепления, изложим вначале для уравнений газовой динамики в декартовых координатах, а затем дадим их обобщение для уравнений Навье—Стокса, в том числе и для уравнений в произвольных криволинейных координатах. Пусть
ди
N
ТО",
™ = -£
3 = 1
дх3
Р = Р(Р,е)
(1)
— система уравнений Эйлера в декартовых координатах. Здесь и, W3• — соответственно векторы искомых функций и потоков в направлении ,. Ограничимся для простоты изложения двумерным случаем. Тогда
(т \
рУзу + 83-р
и
( Р \
ру 1 рУ 2
Vе )
W1
ру 2У + 832р
V V (Е + р)
Е
Р
V2
е
Наряду с консервативной формой запишем исходные уравнения (1) в недивергентном
д£
д
£ ,Р
3=1
р(р,е)
(2)
2
и предельно дивергентном виде
I = -А^.
А = ^, р = р(р,е).
В дальнейшем, для определенности, будем использовать уравнение состояния идеального газа р = (7 — 1)ре. Очевидно, что выбор вектора искомых функций ¥ задает форму матричных операторов В3, что дает возможность варьирования при построении разностных схем. Зададим, например, вектор ¥ в виде ¥ = (р, рг^рг2,р)т. Тогда
В,
(
д33 0
0 0
0
д
д33 0 д1
0
д
дх 3
з
д1
0
л-®
3 дх ^
3 дх д
-— с
- V,
\
3
3
дх, р 3 дх, р дж3 /
здесь с, = 7Р. Данный выбор искомых функций ¥ в (2) позволяет сохранить представление уравнения неразрывности и движения в консервативной форме, что предпочтительно при решении задач газовой динамики и, кроме того, позволит сократить число операций при реализации схемы.
Представим матричный оператор В3 в виде суммы двух операторов (специального расщепления по физическим процессам):
В
в2 + в2
(3)
где
в1
0 0
0
0
0 0
0 д
0 0
0 д
0
3 дж3
3 дх 3 д
-- с
- V,
у 3 дж3 р 3 дж3 р дж3 у
В2
(
0
V
0 0
д
д33 0
0
д
д33 0
0
0
0 0
/
Расщепление операторов В3 на их сумму в форме (3) оптимально, так как приводит к минимальному числу дополнительных членов, возникающих в схемах расщепления при аппроксимации уравнений. Действительно, диссипативный матричный оператор
В1В2
0 0 0
V
0 с
0 0 0
д1д
3 дх3- р дх3 3
0 0 0
д1д
V 3 с
0 0 0
3 дх3 р дх3 3
0
/
0
с
3
0
0
0
1
1
1
1
возникающий в разностных схемах при использовании расщепления, содержит дополнительные члены только в уравнении энергии. Все другие расщепления при сохранении
свойств безусловной устойчивости разностных схем содержат большее число ненулевых членов в диссипативных матрицах В^В2.
Введем в расчетной области, в которой отыскивается численное решение, разностную сетку с равномерными шагами Н, по каждому пространственному направлению х,. Аппроксимируем дифференциальные операторы д/дх, разностными операторами Л, или Л, с порядком 0(Нк) по формулам (см. также [2]) для к = 1 полагаем:
если ь, < 0, то Л,
Л,-,
если ь, > 0, то Л, = Л,— Л, = Л,+; для к = 2 полагаем (при симметричной аппроксимации производных):
Л,- = Л = (Л - ^- + Л,+)/2.
Здесь Л,± = — Г±1)/Н,, Т±/ = — оператор сдвига. С учетом введенных аппроксимаций разностные операторы В,^ аппроксимируют исходные операторы В,:
0 0
0 0
0 0 0
-з з рП
0 0 0
^Л,-— с?Лз —
0
3 3 р«
¿2Лз
/
В
/ 0 0 0
0 0 0
0 0 0
V 0 0 0 0
(4)
с порядком 0(Нк). Для получения безусловно устойчивых схем при аппроксимации конвективных членов и членов с давлением выбраны сопряженные аппроксимации Л, и Л, соответственно, подобно [2]. В дальнейшем при построении разностных схем аппроксимацию первых производных в операторах В,^ будем задавать несимметричной с учетом знака скорости с первым порядком, а в операторе правых частей W« — симметричной со вторым порядком. Для численного решения системы уравнений (1) рассмотрим разностные схемы:
1) схему приближенной факторизации
П(1 + таВ1Л)(1 + таВ]н) ,=1
£«,+ 1 _
т
-(А-1)^«
(5)
или эквивалентную ей схему в дробных шагах:
С
_(А;
(I + ™Вис+1/4 = £«, (I + таВ^С «+1/2 = £ «+1/4, (I + таВ2^)£ «+3/4 = £ «+1/2,
(I + таВ2^)£«+1 = £«+3/4,
£«+1 = £«. + т£
«+1
2
2
2) схему предиктор-корректор:
(I + таВ^):Г+1/8 = Г, (I + таВ?^ п+2/8 =: п+1/8,
(I + таВ2^)Г+3/8 = Г+2/8, (7)
(I + тав2^)г+1/2 = г+3/8,
ип+? — ип
п+1/2
т
Разностная схема (6) аппроксимирует исходные уравнения с порядком 0(т + ^2), а схема (7) — с порядком 0(т2 + Л2) при а = 0.5 + 0(т). Как показывает линейный анализ устойчивости, обе схемы безусловно устойчивы при а > 0.5. Разностная схема (7) консервативна, так как она аппроксимирует систему уравнений (1) на этапе корректора в дивергентной форме. В отличие от нее, разностная схема (4) консервативна лишь при установлении, так как аппроксимирует в дивергентной форме стационарные уравнения = 0, в силу того, что (А-1)"" = 0. Реализация схемы приближенной факторизации (5) или (6) на дробных шагах при нахождении невязок £га+1/4 с учетом аппроксимации операторов В^ из (4) сводится к трехточечным скалярным прогонкам. Подобным образом реализуется и схема (7) на этапе предиктора при вычислении вектора :п+1/8 (/ = 1,..., 4) на дробных шагах. Опишем их реализацию на примере схемы (7). На первом (^ = 1) и третьем (^ = 2) дробных шагах решение разностных уравнений
(I + таВ]^^-1^ = Г+°-1)/4, как следует из вида матриц В*^ сводится к решению системы уравнений:
рП+(2О-1)/8 = рП+О-1)/4,
рп+0-1)/4(^+(2о-1)/8 — ^+°-1)/4) + та^А^^-^8 = 0, рп+0-1)/4(^+(2о-1)/8 — ^+(з-1)/4) + та^2Л3рга+(2з-1)/8 = 0, рп+(2о-1)/8 + та^Л3рп+(2о-1)/8 + тасО Л^п+(2з-1)/8 = рп+О-1)/4.
(8)
Из первого уравнения схемы (8) следует, что значение плотности на новых шагах сохраняется с предыдущего шага. Для нахождения давления на новых по времени слоях
, (п • 1\/о п+(2О-1)/8
п + (2^ — 1)/о исключим компоненту скорости V из последнего уравнения систе-
мы (8). В итоге получим разностные уравнения для давления (для первого дробного шага при ] = 1 и третьего — при ] = 2):
1+та^Ло — т2а27ргаЛо п+(1-1)/4Л^ рп+(2о-1)/8 = рп+(о-1)/4 — та7рпЛ^п+(°'-1)/4.
р
Их решение находится трехточечной скалярной прогонкой, после чего явно вычисляются компоненты скорости г™+(2,7 1)/8 и соответственно импульсы (ргО-)га+(2°-1)/8 по каждому направлению. На четных дробных шагах п + ]/4 (на втором — при ] = 1 и четвертом — при ] = 2) система разностных уравнений
(I + таВ^Р^4 = Г+(2о-1)/8
в скалярном виде примет следующий вид:
рп+,/4 + таК, упрп+3/4 = рп+(2,-1)/8,
(ру1 )п+з/4 + таЛ, уп(ру1)п+з/4 = (рУ1)п+(2з-1)/8, (ру2 )п+з/4 + таЛ, уп(ру2)га+з/4 = (рУ2)га+(2з-1)/8,
(9)
р
,п+3'/4 = р«+(2^-1)/8.
Решение первых трех уравнений из (9) может быть получено скалярными трехточечными прогонками независимо для каждого уравнения, а значение давления на новом слое переносится с предыдущего слоя. Наконец, из последнего уравнения схемы (7) на этапе корректора явно определяются значения всех функций на новом временном шаге. Реализация краевых условий схем на дробных шагах (8),(9) и на полном шаге для схем (6), (7) проводилась подобно [3]. Таким образом, решение системы уравнений газовой динамики по схемам предиктор-корректор сводится к скалярным трехточечным прогонкам. Отметим, что число прогонок по каждому направлению совпадает с числом уравнений в отличие от ранее предложенных схем расщепления [3], что дает экономию по числу арифметических операций на реализацию схемы (при переходе с п-го слоя на п + 1) на 33 % в двумерном случае и 25 % — в трехмерном. В силу симметричной аппроксимации операторов Л,Wj в схемах (5) или (7) в численном решении возникают осцилляции. Для их подавления вводится сглаживающий оператор по каждому направлению второго порядка малости, подобно [8], имеющий следующий вид (например, в направлении 1):
Л^г
Яг+1 — Яг-1 2Н1
— а • sign(v3• )е1-
2 Яг+1 — 2Яг + Яг-1
2Л.1
|Яг+1 — 2Яг + Яг-1| |Яг+1 — Яг| + |Яг — Яг-1|' 0 , если |Яг+1 — Яг| + |Яг — Яг-11
(10)
0.
Отметим, что введение сглаживания по каждому направлению вида (10) приводит к некоторому понижению устойчивости схемы.
Как показывают результаты численных расчетов, применение схемы предиктор-корректор (7) предпочтительнее для решения нестационарных уравнений, а схемы приближенной факторизации (5) или (6) — при решении стационарных задач методом установления.
2. Разностные схемы для уравнений Навье—Стокса
Представим уравнения Навье—Стокса сжимаемого теплопроводного газа в декартовых координатах в векторном виде:
N
^ = w + а w = — V
Г1+ -'
3=1
дх,
(11)
здесь и, Wj — соответственно векторы искомых функций и потоков в направлении
и
р
рУ1 рУ2
Vе У
w,
( ру,
ру1у, + , — а]
рУ2У, + 83-р — а2
2
у, (Е + р) + Я, — X] уг а
\
V
з ^ г=1
2
£
1
3
а,
• • дТ
5!-лV + + —-); = -ко
дж,
дж/
Л = — 2/3^;
к = к0/ср, а к0, , ^ — коэффициент теплопроводности и коэффициенты первой и второй вязкости соответственно, которые будем полагать функциями температуры; ср — коэффициент теплоемкости при постоянном давлении; О — вектор внешних сил. Для замыкания системы уравнений (11) заданы уравнения состояния в следующем виде:
Р = е = е(р,т )•
Как и для уравнений газовой динамики, наряду с дивергентной формой приведем уравнения в недивергентном
д ¥ 2 д£ = — ^ в
(12)
или предельно дивергентном виде
I = —А
А
д и
•
(13)
Использование потоковых переменных для уравнений Навье—Стокса неудобно из-за наличия в них вторых производных от компонент скоростей (элементов тензора вязких напряжений а,) по направлениям ж,, поэтому обычно выбирают примитивные газодинамические переменные. Пусть ¥ = ,р)т — вектор искомых функций. Тогда
д
(
в,
V
0
0 0
' дж,
д
с
0
р дж, 3 дж,
,1 -д
дж
0 0
А 'дж,-
р дж, 1 дж, д
дж,
А
' дж,
0
д
' дж,
' дж,
дд -к-
\
1
дж, дж, р /
где = (Л + + а, = р, с, = 5'7р. Матричные операторы В, содержат газодинамические члены уравнений (конвективные члены и члены с давлением) и часть вязких членов (только повторные производные) по каждому направлению ж,, а вектор Е — оставшиеся члены (смешанные производные в уравнениях движения, диссипативные члены и внешние источники О в уравнении энергии). Аппроксимируем матричные операторы В, разностными операторами В,^ с порядком 0(Л,к)и подобно (3) представим их в форме расщепления
В
V
0 0 0 0
0 0 0 а)Л,-
0 0 0 а?Л,-
0 с1 Л с?Л • V-Л, — Л, К
\
3 р—
( Л,- V 0
—
0 0
0
V-Л, — — Л, Ь1Л, , л р— л ,
0 0
0 0
V—Л,- — —Л,- Ь?Л
р-0
0 0
0 0
V
а
1
0
V
а
с
V
1
1
—
Разностная схема
2 ;п+1 _ ; п
П (I + таБ1л)(1 + таБ^)-т-= — (Л-1)^ + СП) (15)
3=1 т
или эквивалентная ей схема в дробных шагах
еп = — (л-1)^ + сп), (I + тавиеп+1/4 = еп,
(I + тавуеп+1/2 = еп+1/4, ()
(I + тав2^)еп+3/4 = еп+1/2, ( )
(I + таВ2^)£ п+1 = еп+3/4,
;п+1 = ;п + т^п+1
аппроксимирует уравнения Навье—Стокса с порядком О(т + Л2) и в линейном приближении подобно схеме (5) безусловно устойчива при а > 0.5. При установлении схема (15) консервативна. Реализация схемы на дробных шагах сводится к трехточечным скалярным прогонкам подобно схеме (6). Отличие от нее состоит в том, что в матричные операторы В3^ добавлены вторые (повторные) разностные производные, аппроксимируемые также на трехточечном шаблоне.
Аналогично схеме (7) для численного решения уравнений Навье—Стокса может быть построена и разностная схема типа предиктор-корректор:
(I + таВ^Г п+1/8 = р,
(I + таВ^Г п+1/2 = р+1/4,
(I + таВ2^)Гп+3/4 = ;п+1/2, (17)
(I + таВ^)^1 = Р+3/4,
ип+1т— ип = wn+1/2 + сп.
В отличие от схемы (16) она аппроксимирует исходные уравнения (11) с порядком О(т2 + Л2) при а = 0.5 + О(т) и является консервативной при решении как стационарных, так и нестационарных задач. На этапе предиктора на дробных шагах она реализуется скалярными прогонками подобно схеме (16), а на этапе корректора — явно. Для погашения осцилляций в схемы (16) и (17) в векторы Wn и Wn+1/2 вводятся сглаживающие члены вида (10) с порядком О (Л2).
3. Преобразование координат
Разностная схема в новых координатах
При проведении численных расчетов равномерные сетки используются сравнительно редко. Это связано с необходимостью сгущения сеток в областях больших градиентов для повышения точности расчета или в силу других причин. Наиболее распространен метод преобразования координат, позволяющий перевести исходную расчетную область в единичный квадрат и полосу, а необходимые сгущения координат в нужных подобластях задать в преобразовании координат (см., например, [7]).
Введем в физической расчетной односвязной области П, в которой будем искать численное решение, невырожденное преобразование координат
93 = 93
где г = 1, 2. Для него существует обратное преобразование координат
(18)
ж3 Ы.
С учетом преобразования (18) систему уравнений Навье—Стокса (10) в новых переменных 93 вновь можно записать в консервативной форме, в виде (см. [2, 8])
N
ди=w, w=—у ,
3=1
(19)
здесь
и = 1
Vе У
3 =
( д51 \ дх1 дх1
д91 д<?2
12 91
91 91
9192 — 9292
V дж2 дж2 /
— соответственно вектор искомых функций и 3 — якобиан преобразования координат, а
W• = 1 ™3 3
/^3
р^и + 9)р — С1 р^2 ^3 + 92Р — С2
\
V (Е + р)^ + ^ — Д У
— векторы потоков в направлении 93, где соответствующие коэффициенты равны
о о
дТ
и
О А
2
I
1=1
Е93 , с3 = £ 9>т
т=1
93
д-91 Д.
дж3' 3
к
1=1
дг3
дф 3 д9г
Е
1=1 2 2
1
93 °3 = Е
т=1
дг^
Л V + ^Е Ы^ + 917^ , ¿IV V = Е Е^-^3, А = ^ — 2/3^
=11=1 3 д91
Так как сетка в старых переменных после преобразования координат становится неравномерной, то численное определение коэффициентов 93 = д 91^x3 неудобно, поэтому обычно используют коэффициенты ж3 = дж1/д93-, воспользовавшись связью между старыми и новыми координатами: 93 = (—1)!+33ж3, где 3 = 1/(ж1ж2 — ж^!).
Представим систему уравнений Навье—Стокса (17) в не дивергентном виде:
дГ 2
д*
+ Е Б3; = г.
3=1
(20)
Здесь матричные операторы Б3 имеют конвективные члены, члены с давлением и часть вязких членов, содержащих лишь повторные производные по каждому направлению 93,
2
з
а вектор Е — оставшиеся члены уравнений: пусть ¥ = (р,^ , ^2,р)т — вектор искомых функций. Тогда матричные операторы В3 для данного вектора ¥ примут вид
В,
0
0
1 д
с--
3 дЪ
- С2 > 3
С2А 3 %
5
А ¥
- с
3 %
где
с3
■Е
т=1
5 ■ д
Е
т=1
д
5 1
■ Р", ) ' 7
Введем в новых переменных д3 равномерную разностную сетку и аппроксимируем исходные операторы В3 разностными операторами В3-^ с порядком 0(Л,к). Введем подобно (3) расщепление операторов В3^ на сумму операторов В^:
/ 0 0 0 0 \
0 0 0
0 0 0
\0 с)Л3 с2Л3 игаЛ3 - с^/
а]А3 а2л3
В
3^
( Л3 и; 0
0 0
0
0 0
0 0
- с2
_д
3^
3
0
0 0
0 0
(21)
Отметим, что вид матричных операторов В^ в новых переменных (21) подобен виду этих же операторов в декартовых координатах. Тогда для численного решения уравнений Навье—Стокса (19) в новых переменных можно использовать разностные схемы (16) или (17), причем их реализация на дробных шагах с точностью до коэффициентов совпадает с реализацией схем в декартовых координатах.
0
0
0
0
а
0
а
0
1
с
з
Р
2
4. Результаты численного моделирования 4.1. Одномерные тестовые расчеты
Для апробации предложенного численного алгоритма и оценки его эффективности проведены расчеты одномерных и двумерных течений в приближении уравнений газовой динамики и уравнений Навье—Стокса сжимаемого теплопроводного газа. Расчеты подтвердили теоретические оценки по устойчивости схемы и ее достаточной точности при аппроксимации исходных уравнений со вторым порядком по времени (для нестационарных задач) и пространству. В качестве одномерных задач рассматривались две задачи: находилось
у
У
05
0.4
0 9
0.8
0.7
0 6
\
._I_I_I_I_I_I_1_I_,_I_I_I_I_1_I_I_I_I
о
0.25
0.5
О 75
x
О
О 25
0.6
0.75
x
Рис. 1. Распределение плотности в канале
Рис. 2. Распределение скорости в канале
нестационарное решение задачи о распаде произвольного разрыва и методом установления отыскивалось стационарное решение задачи о квазиодномерном течении в канале переменного сечения. Обе задачи имеют точное решение и могут служить тестами для проверки алгоритма. Как показали результаты расчетов, численные решения, полученные по предложенным схемам, с достаточной точностью совпадают с точными решениями. Для примера приведены распределения плотности и скорости в канале переменного сечения после установления (рис. 1 и 2). На входе в канал задавался дозвуковой поток, течение подбиралось таким образом, чтобы внутри канала формировался скачок. Несогласованная аппроксимация производных в предложенных схемах (на дробных шагах — с первым порядком и в правой части — со вторым) приводила к некоторым ограничениям на временной шаг. Наилучшая сходимость к стационарному решению достигалась при шагах сетки, соответствующих числам Куранта 1 < К < 7. Как видно из рис. 1 и 2, ударная волна размазывается на две-три ячейки (сплошная линия — точное решение, кружочки — расчет, все параметры потока отнесены к их значениям на входе в канал ро,Уо,ро). Обе схемы в силу консервативности правильно передают скорости движения фронтов, а введенное сглаживание практически устраняет осцилляции решения, присущие схемам второго порядка аппроксимации.
4.2. Двумерные расчеты. Течение в канале
В следующей серии расчетов приведены некоторые результаты численных расчетов внутренних течений вязкого теплопроводного газа в канале типа "воздухозаборник". Расчетная область представляет собой прямоугольник размером 6.5 (длина) на 2 (высота), внутри которого помещен плоский канал (рис. 3).
Его нижняя стенка, отнесенная от границы расчетной области на расстояние а = 0.5, начинается в точке 0.0, а верхняя снесена по оси абсцисс относительно нижней стенки на 0.5. Высота канала составляет 2.0, длина — 6.0. Решение задачи находилось в приближении уравнений Навье—Стокса вязкого сжимаемого теплопроводного газа (19) в преобразованных координатах (18) по предложенной выше разностной схеме (15), (16) с расщеплением
Рис. 3. Расчетная область
операторов в форме (19). Преобразование координат задавалось в следующем виде:
У1(х0 - 1
51 = (т-) , 52
Ь — а
У2(х1) — У1 (х1)
(22)
что позволяет сгущать узлы расчетной сетки в области пограничного слоя. Здесь у1 (х1), у2(х) — уравнения нижней и верхней стенок канала; показатели степени 9, в — коэффициенты сгущения узлов сетки к границам расчетной области; а — левая граница канала; Ь — правая граница канала.
На входной границе области при х1 = —0.5 задавался невозмущенный сверхзвуковой поток, параллельный оси х1:
на линиях х2 = 0, —0.5 < х1 < 0.0 и х2 = 2, —0.5 < х1 < 0.5 — условия симметрии;
дТ
на стенках канала условия прилипания и тепловой изоляции V, = 0, —— = 0,
дх2
дт$ ^
а на выходе из канала — мягкие условия —— = 0. При проведении расчетов полагали
дхт
т =2.
Для удобства численного интегрирования и анализа результатов исходные уравнения Навье—Стокса (19) были приведены к безразмерному виду, где в качестве безразмерных параметров выбраны следующие величины: высота половины канала Ь, скорость набегающего потока [70, плотность газа р0, температура То. Тогда система уравнений (19) в криволинейной системе координат в безразмерном виде сохраняется в том же виде, а при вязких членах в уравнении движения и в уравнении энергии появляются безразмерные параметры — числа Рейнольдса И,е и Прандтля Рг (см., например, [11]).
Стационарное решение задачи находилось методом установления. Расчеты проводились на равномерных и на неравномерных сетках (при больших числах И,е) при сгущении узлов сетки в области пограничного слоя. В первой серии расчеты проведены при Ы = 2, И,е = 103, 7 = 1.4, Рг = 0.72, и = 0.76 на равномерной сетке, содержащей 521 х 161 узлов. Итерационный шаг задавался равным т ^^ что соответствовало числам Куранта 2 < К < 4. Установление к стационарному решению с точностью (рга+1 — рга)/трга < 10-4 достигалось за ^ 2500 итераций. На рис. 4 представлены распределения скорости нормы вектора скорости, а на рис. 5 — распределение плотности газа в расчетной области. Из рис. 4 хорошо видно, что на стенках канала развивается пограничный слой, а в центре
канала остается сверхзвуковое ядро течения. Наиболее разогретый газ находится около стенок канала, при продвижении газа вдоль границ канала происходит дальнейший прогрев зон, расположенных около пограничных слоев.
В сверхзвуковой части канала газ остается слабо прогретым относительно невозмущенного потока на входе в канал. Небольшое повышение температуры газа возникает в областях прохождения ударных волн. Изолинии плотности в канале представлены на рис. 5. Из рис. 4 и 5 хорошо видно первую Х-образную структуру скачков уплотнения, возникающих на стенках у входа в канал за счет торможения газа на них. Они взаимодействуют и отражаются от пограничного слоя около стенок канала. После взаимодействия с пограничным слоем интенсивность ударной волны падает, и далее в канале образуется достаточно сложная структура взаимодействующих волн. Угол наклона головной ударной волны составляет в к 36°, что соответствует теоретической оценке для данного числа Маха с учетом толщины пограничного слоя (см. [11]).
С целью исследования точности получаемого решения и оценки его сходимости алгоритмов были проведены расчеты на последовательности сгущающихся сеток. Для расчетов была выбрана базовая сетка 521 х 81, которая последовательно дробилась в 2 и в 4 раза. На рис. 6 показаны распределения величин нормальной компоненты скорости в продольном сечении вдоль координаты х2 = 0.25, на рис. 7 — распределения плотности в сечении х2 = 0.125 (решения, полученные на сетках 81, 161 и 321 узел, отмечены на графиках цифрами 1, 2 и 3 соответственно).
На рис. 8 и 9 приведены распределения числа Маха на последовательности сеток в поперечных сечениях XI = 1.125 и XI = 2.25: рис. 8 соответствует течению, частично прошедшему через первый скачок уплотнения, а рис. 9 — течению, прошедшему через первый скачок уплотнения и частично через второй скачок. Из рисунков видно, что решения на
Рис. 6. Распределение компоненты скорости У2
Рис. 7. Распределение плотности
---—.—----
1 75 1,5 г
' г -3
0.76 1
0.5 -0.25 -
Рис. 8. Распределение числа Маха Рис. 9. Распределение числа Маха
грубых сетках имеют меньшую точность, хотя количественные характеристики течения до и после скачков уплотнения совпадают у всех трех решений.
Из приведенных выше результатов можно сделать следующие выводы:
1) предложенная схема второго порядка точности позволяет рассчитывать сложные течения без возникновения у них нефизических осцилляций;
2) разностная схема устойчива в широком диапазоне итерационных параметров т;
3) решения на сетках, содержащих более 161 и 321 узлов в нормальном направлении (т.е. содержащих более 6... 8 узлов в области скачков), достаточно хорошо совпадают. Такие сетки могут быть использованы при проведении расчетов и для получения как качественной картины течения, так и количественных характеристик потока.
4.3. Течения газа в канале с источником вдува газа
С целью проверки работы алгоритма при расчете более сложных течений были рассчитаны течения вязкого газа в канале с источником вдува газа с трансзвуковой скоростью с части его поверхности. Для этого на нижней стенке канала введем источник подачи газа, начинающийся в точке х = 1.0, шириной 0.05 единицы. Данная задача — это упрощенный аналог задачи о подаче топлива в газообразном состоянии с околозвуковой скоростью в канал двигателя. Стационарное решение задачи (19) с соответствующими краевыми условиями находилось при тех же параметрах набегающего потока, что и в предыдущей задаче, а скорость вдува газа полагалась звуковой, т.е. М = 1.0. Ниже приведены результаты расчетов для чисел Ке = 104 и Ке = 103. На рис. 10 и 11 представлены распределение нормы вектора скорости и изолинии давления для числа Ке = 104, на рис. 12 и 13 — для числа Ке = 103. На рисунках хорошо видно, что волновая картина течения существенно изменяется при изменении числа И,е.
Перед источником возмущения появляется ударная волна высокой интенсивности, она взаимодействует с противоположной стенкой и вызывает отрыв пограничного слоя у верхней стенки канала. На рис. 14 и 15 приведены поля скорости около верхней и нижней стенок канала для Ке = 104, на рис 16 и 17 — для Ке = 103. Отметим, что при числе Ке = 103 перед источником вдува газа образуется дополнительно два вихря, отсутствующие для Ке = 104. На верхней стенке канала в окрестности точки взаимодействия ударной волны с пограничным слоем возникает характерная система из двух волн уплотнения и волны разрежения между ними.
Расчеты течений в канале, в том числе и со вдувом газа с части его поверхности, на неравномерных сетках 521 х 161 со сгущением узлов в области пограничного слоя практически совпадают с расчетами на равномерных сетках с числом узлов 521 х 321. Отметим,
Рис. 10. Поле скорости (Ие = 104)
Рис. 11. Изолинии давления (Ие = 104
Рис. 12. Поле скорости (Ие = 103)
Рис. 13. Изолинии давления (Ие = 103)
Рис. 14. Поле скорости, нижняя стенка, (Ие = 104)
Рис. 15. Пограничный слой, верхняя стенка, (Ие = 104)
что введение преобразования координат вида (18) дает увеличение числа членов в уравнениях Навье—Стокса (19) и, как следствие, — увеличение числа арифметических операций на узел сетки. Задача по оптимизации затрат ресурсов ЭВМ на решение приведенных выше задач за счет варьирования параметров сетки, временного шага и выбора различных сгущений сеток в работе не рассматривалась.
Заключение
Для численного решения уравнений газовой динамики и системы уравнений Навье—Стокса вязкого сжимаемого теплопроводного газа в работе предложены разностные схемы с оптимальным расщеплением операторов по физическим процессам и пространственным направлениям. Введен оператор сглаживания для получения монотонного решения. Дано обобщение схем на случай криволинейной системы координат. Получены основные характеристики разностных схем, оценены точность расчета и скорость сходимости.
Приведенные результаты моделирования внутренних течений в канале с возникновением зон ударного перехода, их взаимодействием между собой и со стенами канала, возникновением зон отрыва и других особенностей потока характеризуют эффективность предлагаемого метода и его достаточную точность даже без привлечения специальных алгоритмов сгущения сеток. Простота реализации схемы на дробных шагах, возможность изменения параметров потока и шагов сетки в широком диапазоне (на основании большого запаса устойчивости схемы), экономичность алгоритма делают его эффективным инструментом для исследования проблем аэродинамики.
Список литературы
[1] Флетчер К. Численные методы в динамике жидкостей: В 2 т. М.: Мир, 1991.
[2] Численное решение многомерных задач газовой динамики / С.К. Годунов, А.В. Забродин, М.Я. Иванов и др. М.: Наука, 1989.
[3] Ковеня В.М., Яненко Н.Н. Метод расщепления в задачах газовой динамики. Новосибирск: Наука СО, 1981.
[4] Computer and Fluids, Special Issue, 6th International Symposium on CFD. 1998. Vol. 27, N 5-6.
[5] Computational Fluid Dynamics Journal, Special Issue dedicated to prof. H. Daiguji. 1999. Vol. 8, N 2.
[6] Computational Fluid Dynamics Journal, 4th Asian Workshop on CFD. 2004. Vol. 13, N 2.
[7] Computer and Fluids, Special Issue dedicated to prof. R. Peuret. 2002. Vol. 31, N 4-7.
[8] Ковеня В.М. Разностные методы решения многомерных задач: Курс лекций. Новосибирск: НГУ, 2004.
[9] Ковеня В.М., Лебедев А.С. Численное исследование отрывного течения в ближнем следе. Новосибирск, 1987 (Препр. ИТПМ СО АН СССР. № 14-87). 48 c.
[10] Ковеня В.М., Лебедев А.С. Модификации метода расщепления для построения экономичных разностных схем // Журн. вычисл. математики и мат. физики. 1994. Т. 34, № 6. C. 886-897.
[11] Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1970.
Поступила в редакцию 19 января 2007 г., в переработанном виде — 2 февраля 2007 г.
