Неявный метод численного моделирования пространственных течений вязкого газа Текст научной статьи по специальности «Физика»

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

Том 8, № 1, 2003

НЕЯВНЫЙ МЕТОД ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ПРОСТРАНСТВЕННЫХ ТЕЧЕНИЙ ВЯЗКОГО ГАЗА*

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

The numerical method of solution of full system of Navier-Stokes three-dimensional equations circumscribing viscous gas flows is offered. In the basis of algorithm the final volumes implicit method is fixed. Nonviscous streams are approximated within the TVD-philosophy with use of van Leer limiters. The implicit operator transformation is carried out by means of LU-factorization. The calculations of two-dimensional and three-dimensional problems of viscous gas solid flow are conducted for method testing. The comparison with experimental data and data of other authors's calculations demonstrate good accuracy of offered method for such features flow modelling as boundary layers, separation and association stream.

Введение

За последние 15 лет предложено большое число методов численного решения полной системы уравнений Навье — Стокса, основанных на концепции TVD-монотонизации [1-5]. Эти методы можно условно разделить на два класса по типу аппроксимации невязких потоков: центральноразностные [5] и смещенные против потока [1-4]. Для противопо-точных схем существуют различные способы конструирования потока через грань ячейки. Наряду с традиционным способом, где монотонизации подлежат псевдохарактеристические переменные или сами потоки [4, 5], необходимо отметить MUSCL-подход, применяющийся для противопоточных схем. Здесь поток через грань ячейки вычисляется как функция от вектора примитивных переменных, "восстановленного" на грани по значениям из соседних ячеек с использованием монотонизированной противопоточной интерполяции высокого порядка. Такой подход пременяется в работах [1-3]. Порядок аппроксимации невязких членов в методах [4, 5] - второй, в [2, 3] - третий. Алгоритм Дайгужи [1] обладает четвертым или даже пятым порядком пространственной аппроксимации конвективных членов.

* Работа поддержана Российским фондом фундаментальных исследований, грант 01-01-00799, и грантом конкурса молодых ученых, посвященного 45-летию СО РАН.

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

Все перечисленные выше решатели, за исключением [1], построены на идеологии метода конечных объемов, что автоматически гарантирует консервативность этих алгоритмов. Неявные методы можно дифференцировать также по способу обращения неявного оператора: метод попеременных направлений [2, 4], LU-факторизация [1], релаксация Гаусса — Зейделя [3].

Настоящая работа продолжает цикл работ [6, 7], посвященных моделированию течений вязкой несжимаемой жидкости и распространяет предложенные в них идеи на случай расчета сжимаемых течений.

Характерная особенность сжимаемых течений при сверхзвуковых скоростях — наличие ударных волн — больших градиентов в решении. Поэтому при расчете сжимаемых течений методами повышенного порядка аппроксимации использование монотонизации играет решающую роль. В статье [8] предпринята попытка проанализировать основные противо-поточные и центральноразностные TVD-схемы Хартена — Йи и Чакравати — Ошера с точки зрения их точности и скорости сходимости к стационарному решению. На основании проведенных в этой работе численных экспериментов сделан вывод, что при расчете разрывных решений наилучшим образом сочетает в себе точность и скорость сходимости TVD-схема Чакравати — Ошера второго порядка с "гладкой" функцией-ограничителем Ван-Лира. Поэтому для аппроксимации конвективных членов в настоящей работе берется именно эта комбинация.

Итак, в работе предложен неявный численный метод решения стационарных уравнений Навье — Стокса для расчета ламинарных течений вязкого газа. Алгоритм основан на неявном методе конечных объемов с TVD-аппроксимацией конвективных членов и центральноразностной аппроксимацией вязких членов. При вычислении монотонизирующих поправок используется ограничитель Ван-Лира. Метод обладает вторым порядком аппроксимации по пространству и первым — по времени. Использование в качестве искомых простейших переменных p,u,v,w,T позволяет учесть в неявном операторе якобианы вязких потоков, что повышает устойчивость метода. Для обращения неявного оператора используется метод LU-факторизации, позволяющий свести переход на новый слой по времени к двум обходам расчетной области бегущим счетом.

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

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

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

Вместе с тем в работе проведено сравнение результатов расчетов рассмотренных тестовых задач с данными, полученными другими исследователями с использованием известных алгоритмов CFL3D [2] и LAURA [5]. Показано, что предложенный численный метод

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

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

1. Основные уравнения

Пространственные течения вязкого газа описываются системой уравнений Навье — Сток-са, которая в декартовой системе координат х, у, г имеет вид

dQ + dE + dF + dG =0 dt dx dy dz

(1)

где

Q

f p\

pu pv pw

W

E

pu

pu2 + p - ßTxx puv - ßTXy puw - ßTxz \(e + p)u - pbi - kTxJ

(

F

pv

\

puv - ßTxy

pv2 + p - ßTyy

pvw - ßTyz

\(e + p)v - ßb2 - kTyJ

G

pw

\

puw

pTxz pTyz

pvw pw2 + p - pTz \(e + p)w - ßb3 - kTzJ

bl = Txxu + Txy v + Txz w b2 = Txy u + Tyy v + Tyz w,

yy

yz

b3 = Txz u + Tyz v + Tzz w

Здесь р - плотность; и,у,,ш - компоненты скорости; р - давление; е - полная энергия на единицу объема; тхх,... ,тхх - компоненты тензора вязких напряжений. Предполагается, что газ - идеальный с показателем адиабаты 7 = 1.4, а коэффициент вязкости р зависит только от абсолютной температуры Т и вычисляется по формуле Саттерленда

р = pr

3/2

Tr + Ts T + Ts

где Tr = 273K, Ts = 120K, pr = p(Tr). Коэффициент теплопроводности k рассчитывается по формуле k = pcp/Pr, где cp - удельная теплоемкость при постоянном давлении, а Pr = const — число Прандтля.

Система уравнений Навье — Стокса обезразмеривается на характерные масштабы длины, скорости, времени, плотности, давления, температуры и вязкости L*,U*,t*, р*,p*,T*, p*. Для сохранения вида конвективной части исходных размерных уравнений полагалось

t*

L*

W*

p*

pÜ.

Характерная температура принималась равной Т = и/су. При этом уравнение состояния в безразмерных переменных приобретает вид

Р = (7 - 1)р'Т'.

Число Рейнольдса задачи определяется из соотношения И,е = Р—-—-. Тогда для пол-

V*

ной постановки задачи и нахождения ее решения в безразмерных переменных достаточно задать число Маха, число Рейнольдса, температурный фактор на теле (отношение ^ где Tw - температура на стенке обтекаемого тела, а To - температура торможения) и размерную температуру набегающего потока на бесконечности TCX). Последний параметр используется только для нахождения безразмерных ^ и ^ в формуле Саттерленда.

2. Численный метод

В основу численного алгоритма положен неявный метод конечных объемов. Система (1) в безразмерных переменных записывается в виде интегральных законов сохранения для произвольного объема V (штрихи опущены):

А

И

J Q ¿V + у & ■ пйБ = 0,

V дУ

(2)

где п - единичный вектор внешней нормали к поверхности дV объема V; & - матрица, составленная из столбцов векторов потоков & = (Е|Ё|0). Она может быть разложена в сумму невязкой и вязкой составляющих

& = &™ + &ь

(3)

/ ри ру рш

ри2 + р риу риш

&т риу риш \(б + р)и ру2 + р руш (е + р)у руш рш2 + р (е + р)ш/

/ 0 0 0

1 ^Тхх ^Тху

рТху Р>Туу

+ р ^Гх 'У ^2 + ^ ^т.у рЬз + р

В расчетной области вводится регулярная сетка с ячейками Vi,j,k в виде криволинейных параллелепипедов. Усредненное по объему Vi,j,к значение вектора неизвестных переменных Q относится к центру ячейки. Тогда для каждой ячейки можно записать следующую неявную аппроксимацию интегральных законов сохранения (2) :

Qij¡к Q

п

ук

V

Vк

К-Н^1-

В правой части (4) собраны потоки через грани ячейки (г,], к), взятые на (п + 1)-м слое по времени:

КНБП+к1

-(Ё+2 - Ё^2 + Ё^1 - Ёj_2 + Ёк+2 - Ёк-2)

п+1

2 "2 + 2 ■> 2 1 2 2 '

Для получения стационарного решения итерации по времени продолжаются до установления.

2.1. Вычисление разностных потоков

Разностные потоки ЁП+1, , аппроксимируют входящие в уравнения (2) интег-

г+ 2 ]+ 2 2

ралы по граням ячейки от потока К. В соответствии с расщеплением (3) представляются в виде сумм невязких и вязких разностных потоков: ]]т+1 = ]]т+ 1 + ]]т+1, т = г,], к.

Для вычисления ]]т+1 применяется противопоточная ТУБ-схема Чакравати — Ошера

[9] второго порядка аппроксимации с ограничителями Ван-Лира. При этом используется расщепление матрицы Якоби вектора потока по знакам ее собственных значений:

]т+1/2 = 2 [(Ит+1 + ^^+1/2 - ^(Л)т+1/2Ат+1/2д] - Лт+1/2, (4)

А = (д(К-■ я)\ =(_ _ Т)

Лт+1/2 =1 -- I = (И ■ Б ■ Т)т+1/2,

V д О /т+1/2 Б = Diag(Л1,Л2,Л3,Л4,Л5), Л1'2'3 = и, Л4'5 = и ± с|Б|.

Здесь и = пБх + уБу + /шБх; п, - компоненты скорости; с - скорость звука на рассматриваемой грани ячейки; Я = (Бх, Бу, Бх)т+1/2 - вектор нормали к грани т +1/2, равный по модулю площади этой грани и направленный в сторону увеличения индекса т; Ит+1/2, Ьт+1/2 - матрицы правых и левых собственных векторов матрицы Лт+1/2, такие что ИТ = I (единичная матрица).

Все необходимые газодинамические параметры на грани ячейки определяются методом усреднения по Роу [10] через значения вектора О в прилегающих к грани ячейках. Функция энтропийной коррекции имеет вид

И, и > ев,

^(Л) = (И- Ь), #*) = < *2 + е2 ,, е8 = е- |Я|. (5)

2 ев

|г| < ев,

Исследование влияния энтропийной коррекции на точность получаемого решения показало, что ненулевая величина е крайне нежелательно сказывается на точности воспроизведения теплового потока и коэффициента трения на поверхности тела при больших числах Рейнольдса. Поэтому при вычислении потоков в нормальном к твердой стенке направлении полагалось е = 0. Необходимость использования е > 0 в касательных направлениях обусловлена следующим. При расчете обтекания затупленных тел касательная к поверхности составляющая скорости в окрестности линии торможения в лобовой части обтекаемого тела стремится к нулю. Если е = 0, то это приводит к занулению собственных значений Л1, Л2, Л3 матрицы ^(Л)т+1/2 в (4). Таким образом, исчезает диссипация базовой схемы первого порядка (соответствующая потоку (4) без члена Л), что приводит к возникновению неустойчивости в окрестности линии торможения. Поэтому при вычислении потоков в касательном к твердой поверхности направлении параметр энтропийной коррекции е принимался равным 0.2.

В выражении (4) Л содержит антидиффузионные члены с введенными в них функциями-ограничителями, обеспечивающие второй порядок аппроксимации по пространству:

Лт+1/2 = 1(Ит+3/2 ' ¿Гт+3/2 - Ит-1/2 ' ¿Г,^-1/2^

Здесь компоненты векторов сгт+3/2 и ¿гт_1/2 вычисляются через

¿^+1/2 = ± |0|)т+1/2Ьт+1/2((т+1 -с использованием функции-ограничителя Ван-Лира ф(х,у) = ^ + I у1 .

¿г

х + у

_ ^т+3/2^ т+1/2 1 I ^т+3/2^ т+1/2 1

ат+3/2ат+1/2 + |ат+3/2СТт+1/2 1

т+3/2 __£ + __£ '

ат+3/2 + ат+1/2

+ 1а+г а+£ I

~ _ "т_1/2^т+1/2 + |"т_1/2^т+1/2 1

¿т_1/2 = +-+€ •

ит_1/2 + °т+1/2

Применение таких ограничителей вместо традиционных минмодульных позволяет существенно повысить скорость сходимости метода установления к стационарному решению [8].

Вычисление вязких составляющих потоков основано на формуле Е= • Я. Аппроксимацию компонент тензора вязких напряжений т и градиента температуры Т, входящих в определение рассмотрим на примере тхх:

4 ди 2 / Зу дт хх 3 дх 3 \ ду дг

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

д д дп д д( д

дх дх д£ дх дп дх д( Входящие в последнее выражение производные по аппроксимируются центральными

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

I

(О! ЗП = _Д£_(я я я )1

V дх'ду'дг Л+1/2 = К+1/2 (Ях'Яу 'Я* )г+1/2'

дп дп дп\ Дп (я п

(ях, яу' яг )г+1/2'

дх ду дг У ¿+1/2 ^+1/2

д(д(д(У = (я я я )К

дх'ду'дг) ¿+1/2 = V+1/2 (Ях ' Яу ' Яг )г+1/2 '

где V+1/2 = 0.5(¥ + ¥¿+1), а векторы (Ях , Яу, Яг)^+1/2 и (Ях , Яу, Яг)К+1/2, направленные в сторону увеличения индексов ] и к, вычисляются усреднением нормалей к соответствующим граням из ячеек с индексами г и г + 1. Аналогичным образом находятся производные на гранях ] + 1/2 и к + 1/2.

2.2. Реализация алгоритма

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

Е ™ = о 1(( ?) + о 2(( п) + С 3(( с).

Система нелинейных разностных уравнений (4) линеаризуется методом Ньютона относи-

П П п.

5 , , Чс •

тельно Qn, Qn Qn Qn

Vj G (д(RHS)V i (Q„+1

AG )(Qj -

Qjfc) dQ a J

^d(RHS^ra(Qa+1 - Qa) = RHsn?fc. (6)

Здесь G = dQ/dQ. Для упрощения неявного оператора при линеаризации вязких потоков в левой части оставляются лишь члены, соответствующие повторным производным вдоль каждого из направлений i,j, k [7]. Иными словами, полагается, что

dG 1

(Fvis )n+1 _ (jpDis )n | I UKJ I _ (Qn+1 — Qn) I /OiVA 4-\2\

Обозначим

,д Q d,

■ (Q n+1 - Q П) + o((At)2),

/S vis \n+1 /S vis \n I (F j+1/2) = (F j+1/2) +

/"S vis \n+1 _ /"S vis \n I

(F fc+1/2,) = (Ffc+1/2) +

д CS 2

д Q n

д G3 д Q с

i+1/2

j+1/2

fc+1/2

(Q П+1 - Q П) + o((At)2),

(Q П+1 - Q n) + O((At)2).

в

i±1/2

в

j±1/2

в

fc±1/2

д G1

30 ^

.^V j±1/2 дG3 \

д Q Jk±1/2

Тогда линеаризованная система (6) может быть представлена в виде (О + Сг_1/2Тг-1 + С;/_1/2Т;/_1 + ^-1/2^-1 +

(7)

+ Ci+1/2Ti+1 + Cj+1/2Tj+1 + Cfc+1/2Tfc+1)A^ij.+1 где T оператор сдвига, Ti+1 ^¿jfc — ^i+1,j,fc;

A^ij+fc = Qn+ - Qjfc;

RHSn/-fc,

IS

Bi-1/2 - A++_1/2Gi-1; ij-1/2 - A+-1/2Gj-1;

Ci-1/2 Cj-1/2 = Bj-Cfc— 1 /2 = Bfc—1/2 - A+-1/2Gfc-1;

^ijfc

Ci+1/2 Cj+1/2 Cfc+1 /2

B,

в

B

+ 1/2 + A—+1 /2 Gi+1; j+1/2 + A-+1/2Gj+1; fc+1/2 + A-+1 /2Gk+1;

Gnik 1 ( A+1 /2 Ai+1 /2 1 A+1 /2 A- + 1 /2 1 A+_1 /2 Al

Gn-G-

A^ijk ^ ^i-1/2 л,+ 1/2 "Г" -г\?-1/2 л^+1/2 ^ "^4-1/2 лк:+1/2у Gijfc - (Bi-1/2 + Bi+1/2 + Bj-1/2 + Bj+1/2 + Bk-1/2 + Bfc+1/2) .

n

n

ЬИ-факторизация. Неявный оператор в левой части (7) приближенно факторизует-ся на два оператора, один из которых содержит все операторы сдвига на узел назад, а другой — все операторы сдвига на узел вперед [6]:

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

2.3. Граничные условия

Для реализации граничных условий за каждой границей вводятся два ряда фиктивных ячеек, значения искомых переменных в которых при I = £п заполняются в соответствии с типом граничных условий. Такой подход позволяет использовать общие процедуры вычисления потоков для всех внутренних, в том числе приграничных, ячеек расчетной области. Выделяются типы границ: входная, выходная, твердая стенка, удаленная граница, плоскость или ось симметрии. Наибольший интерес представляют условия на твердой стенке и удаленной границе, на которых мы и остановимся подробно.

Условие прилипания V = 0 на твердой стенке реализовывалось следующим образом. Скорость в фиктивных ячейках задавалась по принципу антисимметрии так, чтобы среднее арифметическое скоростей из внутренних приграничных ячеек и фиктивных ячеек было равно нулю (рис. 1):

Температура в фиктивных ячейках либо задавалась равной Т^ (в случае обтекания тел с фиксированной температурой поверхности), либо определялась из условия адиабатичнос-ти стенки дТ/дп = 0 в случае обтекания теплоизолированных тел.

х (13 + Сг+1/2Тг+1 + Cj+l/2Tj+1 + Ск+1/2Тк+1)Д^П+ = КН8П

V2 = ^э, V! =

Рис. 1. Определение скоростей в фиктивных ячейках. пС<?

Точное выражение для градиента давления на стенке можно получить из уравнений сохранения импульса с учетом v|r = 0:

öp ön

n

( d|UT5

Re

öx öx

d^Txz \

+ +

\ öx

xz

öz

d^Tyz

öz

d^Tzz

dz /

(8)

В криволинейной системе координат соотношение (8) выглядит довольно громоздко [11]. Как показывают численные эксперименты, использование условия (8) приводит к решениям, в которых на стенке др/дп ~ 0. Поэтому давление в фиктивных ячейках определялось из более простого условия др/дп = 0. При больших числах Рейнольдса (И,е > 1000) вблизи тела справедливо приближение тонкого пограничного слоя, давление в котором не меняется по нормали к поверхности, и условие др/дп|г = 0 имеет строгое обоснование.

В настоящей работе для реализации условий на удаленной границе использован подход, предложенный Томасом и Саласом в [12] для двумерных невязких течений. При переносе этого подхода на вязкие течения предполагается, что влияние вязкости в окрестности удаленных границ несущественно. Пусть п - единичный вектор нормали к границе, направленный внутрь расчетной области. Рассмотрим уравнения газовой динамики (уравнения Эйлера) в проекции на эту нормаль. Собственными числами матрицы Якоби для полученной таким образом одномерной системы будут Л1 = пп, Л2 = пп — с и Л3 = пп + с, где пп = V ■ п, а с - скорость звука. Предположим, что течение в окрестности границы изэнт-

ропическое. Тогда существуют инварианты Римана г2 = пп--- и г3 = пп +

причем

Y - 1

Y - 1

(х

г2 сохраняется вдоль кривой — = Л2,

(И

(х

г3 сохраняется вдоль кривой — = Л3.

(И

Далее, согласно [12], будем предполагать, что энтропия в = р/р7 и касательная составляющая вектора скорости vr = V — V ■ п сохраняются вдоль характеристики (х/(£ = Л1 . Тогда можно сформулировать следующий принцип заполнения фиктивных ячеек на удаленной границе.

Если Л1 > 0 на границе, то г1 и вектор vr в фиктивных ячейках берутся из параметров набегающего потока, иначе г1 и vr экстраполируются изнутри области.

Если Лг > 0, то гг в фиктивных ячейках берется из параметров набегающего потока, иначе гг экстраполируется изнутри области (г = 2, 3).

Таким образом в каждой фиктивной ячейке мы получаем значения комплексов

p/pY

2c

vT

Y - 1

Un +

2c

Y - 1'

1

по которым затем восстанавливаем значения плотности, скорости и температуры.

3. Результаты расчетов

Численный метод реализован в виде пакета программ на языке FORTRAN. Расчет одной итерации для одной ячейки сетки занимает 1.6 • 10-4 с рабочего времени процессора Celeron 533 МГц. Затраты оперативной памяти составляют 30 Мбайт на каждые 100 тысяч ячеек сетки. Для получения установившегося решения требуется провести от одной до десяти тысяч итераций по времени в зависимости от сложности моделируемого течения.

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

1) обтекание цилиндра слабосжимаемым потоком;

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

3) гиперзвуковое течение в плоском угле сжатия;

4) обтекание затупленного конуса под углом атаки 20°.

Основные параметры тестовых задач сведены в таблице. Число Прандтля Pr, если не указано противное, полагалось равным 0.72.

Основные параметры тестовых задач

№ теста Геометрия Re T», K Tw /To Размерность

1 Цилиндр 0.1 40 - адиаб. 2D

2 Ударная волна + пластина 2.0 2.96 ■ 105 и 167 адиаб. 2D

3 Угол сжатия 14.1 1.04 ■ 105 88.9 0.082 2D, 3D

4 Затупленный конус 10.6 3.937 ■ 106 47.3 0.27 3D

Во всех тестах формулы для аэродинамических коэффициентов cp,cf ,Ch,qw приведены в размерных переменных.

3.1. Обтекание цилиндра

Для демонстрации работоспособности алгоритма при малых числах Маха проведен расчет обтекания теплоизолированного цилиндра невозмущенным потоком с = 0.1 Re = 40, Pr = 0.7. Число Рейнольдса посчитано по диаметру цилиндра. В этой задаче вязкость вычислялась по степенной формуле ^ = ß*(T/Т*)ш, где ш = 0.76, поэтому для проведения расчета в безразмерных переменных задание температуры не требуется.

Для сравнения используются данные численных расчетов этой задачи в несжимаемой постановке [13] (М = 0), а также данные эксперимента по визуализации рециркуляционной зоны, проведенного в несжимаемой жидкости [14].

При численном решении задачи внешняя граница области была отнесена на расстояние 20 радиусов от центра цилиндра. Для реализации краевых условий на удаленной границе использовались одномерные инварианты Римана. Ниже приведены результаты расчетов на сетке 80 х 120 ячеек, которые уже не меняются при сгущении сетки (grid independent). Шаг по времени т полагался равным 10.0. На рис. 2 представлено распределение коэффициента давления cp = 2(p — p^)/р^п200 на поверхности цилиндра, полученное описанным выше методом, в сравнении с расчетом из работы [13], где решались несжимаемые уравнения Навье — Стокса. Наблюдается хорошее совпадение.

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

Рис. 2. Коэффициент давления на поверхнос- Рис. 3. Обтекание цилиндра: история сходи-ти цилиндра: ср = 2(р — р^. мости. Шаг по времени т = 10.0.

На рис. 3 представлена история сходимости к стационарному решению для двух значений числа Маха набегающего потока = 0.1 и = 0.01). Здесь по оси ординат отложена невязка решения стационарных уравнений в норме С:

Егг(п) = шах ( тах

\RHSj [т]

При = 0.1 сходимость до 10_12 достигается примерно за 6000 шагов по времени. При уменьшении характерного числа Маха, как и следовало ожидать в связи с жесткостью исходной системы уравнений, сходимость существенно замедляется. Однако и в этом случае достаточно провести 15 000 итераций по времени, чтобы получить установившееся решение с невязкой 10_8. Этот результат доказывает, что снижение точности метода при М ^ 0, отмеченное выше, целиком обусловлено способом пространственной аппроксимации уравнений и не связано с ухудшением сходимости.

Картина течения в окрестности цилиндра, рассчитанная при = 0.1, показана на рис. 4 вместе с фотографией эксперимента, проведенного в воде при Ие = 41 [14]. Численно найденные точка отрыва потока, длина и форма рециркуляционной зоны с хорошей точностью совпадают с экспериментальными данными.

-1 0 1 2 3 х

Рис. 4. Обтекание цилиндра: изолинии давления и линии тока.

3.2. Взаимодействие ударной волны с погранслоем на пластине

Косой скачок уплотнения, взаимодействующий с пограничным слоем, приводит к отрыву последнего и формированию в месте взаимодействия зоны возвратного течения. Угол наклона косого скачка уплотнения 9 = 32.60. Расчетная область представляет собой прямоугольник размера 1.6Ь х 1.1 Ь. Здесь Ь = 0.049 м — расстояние по оси х от переднего края пластины до точки падения скачка уплотнения, взятое в качестве характерной длины для определения числа Рейнольдса Ие = 2.96 • 105. Расчет проводился на сетке 90 х 168, которая сгущалась к поверхности пластины. На пластине ставилось условие адиабатичности дТ/дп = 0. Значение Тж для этой задачи бралось равным 167 К, что соответствует температуре торможения потока Т0 = 300 К. Было показано, что решение, полученное на сетке 90 х 168, является предельным и не меняется при дальнейшем сгущении сетки.

Результаты численных расчетов сравнивались с экспериментальными данными из работы Хаккинена [15] и расчетами Томаса — Уолтерса [3]. Распределение давления и коэффициента трения на поверхности пластины представлено на рис. 5. Коэффициент трения вычислялся по формуле

2тп, Ч =-- •

Рж и

где ^ = р

дУг дп

. Видно, что давление с хорошей точностью совпадает как с расчетом

^ -___

по схеме Томаса — Уолтерса, так и с экспериментальными данными. Для предложенного метода можно отметить увеличение длины отрывной рециркуляционной зоны (рис. 5, б, область, где 04 < 0) по сравнению с данными Томаса — Уолтерса и экспериментом, а также заниженную величину 04 в задней части рециркуляционной зоны (х/Ь и 1.1). Черными кружками на рис. 5, б обозначены точки, в которых экспериментально зафиксировано отрывное течение, но значения коэффициента трения 04 получены не были. Здесь очевидно общее для двух численных методов занижение 04 за отраженным скачком (х/Ь > 1.2). В целом результаты расчета удовлетворительно согласуются с экспериментом.

Рис. 5. Распределение давления (а) и коэффициента трения (б) на пластине.

3.3. Течение в угле сжатия

Расчетная область и структура течения в угле сжатия в = 24° схематично представлены на рис. 6. Угол сжатия состоит из горизонтальной пластины и следующего за ней клина. В случае достаточно больших углов клина при обтекании этой модели образуется зона возвратного течения в окрестности вершины. Входные данные для этой задачи (см. таблицу) соответствуют эксперименту Холдена и Мозеллы [16], с которым проводим сравнение.

В большинстве работ, посвященных численному исследованию этой задачи, она решается в двумерной постановке. Однако в [17] показано, что конечная ширина модели, использованной в эксперименте, играет существенную роль. Действительно, ширина установки в эксперименте й = 0.6096 м сравнима с длиной горизонтальной части хс = 0.439 м. Поэтому в настоящей работе наряду с двумерными расчетами на сетке 100 х 100 представлены также результаты расчета этой задачи в трехмерной постановке, проведенные на сетке 100 х 100 х 25. Расчетная область для моделирования трехмерного обтекания показана на рис. 7 штриховой линией. Здесь твердая поверхность угла сжатия ограничена

Ударная волна от передней кромки пластины

Рис. 6. Структура течения в угле сжатия.

Р

А

Рис. 7. Геометрия расчетной области, изолинии давления и линии тока на поверхности угла сжатия для трехмерного расчета.

Рис. 8. Коэффициент давления ср = 2р/ржРис. 9. Коэффициент теплопередачи с^ для

для угла сжатия 0 = 24°.

угла сжатия 0 = 24°; ch = qw/ рж(И^ — Hw).

сплошной линией. Разбиение области в поперечном направлении выбиралось так, чтобы на поверхность модели между плоскостью симметрии и боковой стороной (соответствующей точке F) приходились 19 ячеек сетки и 8 ячеек — на область течения между боковой стороной модели и внешним потоком (отрезок FA). Сетка сгущалась к твердой поверхности в нормальном направлении и к переднему краю пластины и вершине угла в направлении оси x. Для реализации условий на границах с внешним потоком, соответствующих граням ABC DE, AFGE, EGHD, использовались одномерные инварианты Римана.

На рис. 8 показано распределение коэффициента давления cp = 2р/ржи^ вдоль поверхности, полученное экспериментально, настоящим методом и с помощью программы CFL3D [17]. Видно, что при решении задачи в двумерном приближении точка отрыва погранслоя xs, которой соответствует скачок давления в окрестности x/xc = 0.5, для обоих численных методов предсказывается неудовлетворительно. Максимум в распределении коэффициента давления сильно смещен вправо. На рис. 9 представлено соответствующее распределение коэффициента теплопередачи ch = qw/ржиж (Нж — Hw). Здесь также наблюдается существенное отличие 2D-расчета от экспериментальных данных.

Распределение cp и ch при решении задачи в трехмерной постановке вычислялось в сечении z = const, проходящем по центру ячеек, примыкающих к плоскости симметрии. Результаты трехмерных расчетов по настоящему методу и CFL3D гораздо лучше совпадают с экспериментальными данными, нежели данные расчетов в двумерной постановке (см. рис. 8 и 9). Таким образом, еще раз подтверждается вывод работы [17] о том, что для проведения корректного сравнения с экспериментом Холдена — Мозеллы необходимо учитывать пространственную структуру течения в угле сжатия. При адекватной постановке задачи предложенный численный метод, равно как и CFL3D, c хорошей точностью передает особенности течения.

Изолинии давления и линии тока на поверхности угла сжатия приведены на рис. 7. Видна тенденция к сужению рециркуляционной зоны с приближением к краю пластины.

3.4. Затупленный конус

Один из наиболее важных показателей численных методов решения уравнений Навье — Стокса — это точность воспроизведения градиентных характеристик в пограничном слое

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

Затупленный по сфере конус с углом полураствора 15° расположен во внешнем потоке под углом атаки а = 20°. Радиус затупления Rn = 0.028 м, радиус основания конуса Rb = 0.1524 м. Параметры набегающего потока и температурный фактор приведены в таблице. Характерный линейный размер для определения числа Рейнольдса — 1 м. Опубликованы данные эксперимента Клэри по определению теплового потока на поверхности конуса [18] и результаты расчета этой задачи с помощью программы LAURA [19].

В настоящей работе расчет течения около конуса выполнен на сетке, содержащей 50 узлов в продольном направлении, 60 узлов между поверхностью конуса и внешней границей, 61 узел в окружном направлении. Сетка сгущалась к поверхности конуса по экспоненциальному закону. Высота ближайшей к поверхности конуса ячейки составляла 0.03 • 10

-4

м

-4

в лобовой части и 0.85 • 10" ' м у основания на подветренной стороне конуса.

Распределение теплового потока на поверхности конуса в сравнении с данными эксперимента и расчетов по LAURA представлено на рис. 10 (вдоль образующих конуса) и 11 (в окружном направлении при фиксированных x). Здесь ф = 0 соответствует образующей конуса, лежащей в плоскости симметрии на наветренной стороне тела, ф = 180° — на подветренной. На графиках размерный тепловой поток qw определяемый как

qw=fc(£) w •

где k — коэффициент теплопроводности, отнесен к qref = 2.15 • 105 МВт. Кривые qw, рассчитанные по настоящему методу, практически совпадают с полученными по программе LAURA в передней части конуса. Особенно хорошо это видно по распределению теплового потока в окружном направлении (рис. 11) при x/Rn = 2.23, где кривые, соответствующие настоящему методу и алгоритму LAURA, накладываются друг на друга. Занижение теплового потока в задней части конуса обусловлено скорее всего недостаточным сгущением сетки к поверхности конуса (сетка для расчета [19] имела 80 узлов в нормальном направлении). В целом же оба численных метода хорошо согласуются с экспериментом.

Рис. 10. Распределение теплового потока вдоль образующих конуса.

Рис. 11. Распределение теплового потока на поверхности конуса в окружном направлении.

плоскость симметрии

присоединение потока в поперечном сечении • расчетной области

линии тока

Рис. 12. Линии тока на поверхности конуса.

Рис. 13. Линии тока на поверхности и в поперечном сечении у основания конуса. (Вид со стороны затупления вдоль оси симметрии конуса.)

На рис. 12 приведены линии тока на подветренной поверхности тела. Сравнение результатов расчетов, выполненных по предложенному алгоритму, с данными, полученными с использованием алгоритма LAURA, позволяет сделать вывод: внешняя граница зоны поперечного возвратного течения для обоих алгоритмов совпадает, однако в настоящих расчетах у основания конуса появляется вторичная рециркуляционная зона, хорошо заметная на рис. 13, в верхней части которого изображены линии тока в поперечном сечении расчетной области.

Заключение

В работе предложен неявный метод конечных объемов для нахождения стационарного решения полной системы трехмерных уравнений Навье — Стокса сжимаемой жидкости. Численный метод имеет второй порядок аппроксимации по пространству. Стационарное решение находится путем установления. Сравнение численных расчетов с экспериментальными данными и результатами расчетов других исследователей позволяет сделать вывод о том, что предложенный метод обладает достаточной точностью расчета течений вязкого газа в диапазоне характерных параметров: 0.1 < M < 15, 40 < Re < 5 • 106. Алгоритм обладает высокой достоверностью определения точек отрыва и присоединения пограничного слоя, правильно передает величину теплового потока при числах Рейнольдса порядка 106. При расчете слабосжимаемого обтекания цилиндра обнаружено, что уменьшение характерного числа Маха течения до значений ниже 0.1 влечет заметное ухудшение точности метода.

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

В дальнейшем планируется провести модификацию метода для проведения расчетов с малыми (M < 0.01) числами Маха и распространить его на случай расчета турбулентных течений.

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

[1] Yuan X., Daiguji H. A specially combined lower-upper factored implicit scheme for three-dimensional compressible Navier-Stokes equations // Computers & Fluids. 2001. Vol. 30. P. 339-363.

[2] Vatsa V. N., Thomas J. L., Wedan B. W. Navier — Stokes computations of a prolate spheroid at angle of attack // J. od Aircraft. 1989. Vol. 26, No. 11. P. 986-993.

[3] Thomas J. L., Walters R. W. Upwind relaxation algorithms for the Navier — Stokes equations // AIAA J. 1987. Vol. 25, No. 4. P. 527-534.

[4] Wada Y., Kubota H. Numerical simulation of Re-Entry flow around the space shuttle with finite-rate chemistry // J. of Aircraft. 1992. Vol. 29, No. 6. P. 1049-1056.

[5] GNOFFO P. A. Upwind-Biased, Point-Implicit Relaxation Strategies for Viscous Compressible Perfect-Gas Flows. NASA TP-2953, Feb. 1990.

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

[7] Черный С. Г., Шашкин П.А., Грязин Ю.А. Численное моделирование пространственных турбулентных течений несжимаемой жидкости на основе k — е-моделей // Вычисл. технологии. 1999. Т. 4, №2. С. 74-94.

[8] Чирков Д.В., Черный С.Г. Сравнение точности и сходимости некоторых TVD-схем // Вычисл. технологии. 1999. Т. 5, №5. С. 87-108.

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

[10] Roe P. L. Approximate Riemann solvers, parameter vectors and difference schemes // J. Comput. Physics. 1981 Vol. 43. P. 337-372.

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

[12] Thomas J. L., Salas M. D. Far field boundary conditions for transonic lifting solutions to the euler equations // AIAA J. 1986. Vol. 24, No. 7. P. 1074-1080.

[13] Квак Д., Чэнг Дж. Л. К., Шенкс С. П., Чакраварти С. Р. Метод решения уравнений Навье — Стокса для трехмерных течений несжимаемой жидкости с использованием простейших переменных // Аэрокосм. техника. 1987. №2. С. 144-153.

14] ВАН-ДАЙК M. Альбом течений жидкости и газа. М: Мир, 1986.

[15] Hakkinen R. J., Greber I., Trilling L., Abarbanel S. S. The Interaction of an Oblique Shock Wave with a Laminar Boundary Layer. NASA Memo 2-18-59W, March 1959.

[16] Holden M. S., Moselle J. R. Theoretical and Experimental Studies of the Shock Wave-Boundary Layer Interaction on Compression Surfaces in Hypersonic Flow. ARL 70-0002, Aerospace Research Laboratories, Wright-Patterson AFB, OH, Jan. 1970.

[17] Rudy D. H., Thomas J. L., Kumar A., Gnoffo P. A., Chakravarthy S. R. Copmutation of laminar hypersonic compression-corner flows // AIAA J. 1991. Vol. 29, No. 7. P. 1108-1113.

[18] Cleary J. W. Effects of Angle of Attack and Bluntness on Laminar Heating-Rate Distributions of a 15° Cone at a Mach Number of 10.6. NASA TN D-5450, Oct. 1969.

[19] Greene F. A. Application of the multigrid solution technique to hypersonic entry vehicles // J. of Spacecraft and Rockets. 1994. Vol. 31, No. 5. P. 744-750.

Поступила в редакцию 28 октября 2002 г.