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

МОДЕЛИРОВАНИЕ ПРОЦЕССОВ

УДК 533.6.011.8

С. Т. Суржиков

МЕТОД РАСЧЕТА СВЕРХЗВУКОВОГО ОБТЕКАНИЯ СФЕРЫ НА ОСНОВЕ АШМ КОНЕЧНО-РАЗНОСТНЫХ СХЕМ

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

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

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

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

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

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

Представленный метод основан на использовании для численного интегрирования уравнений Навье-Стокса ЛиБМ [6] конечно-разностных схем и на применении полностью неявных конечно-разностных схем для численного интегрирования уравнения сохранения энергии, которое формулируется в виде уравнения Фурье-Кирхгофа, т. е. относительно температуры.

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

Постановка задачи. Расчет динамики вязкого теплопроводного газа проводится с использованием уравнений авье-Стокса в двумерной осесимметричной постановке. Расчетная область вкл чает в себя поле течения в невозмущенном потоке, за фронтом ударной волны, в следе за обтекаемым телом. Течение предполагается ламинарным во всей области.

Исходная система уравнений формулируется в следующем виде:

^ + div(pV) = 0;

(1)

t + div (-V) = -1 - 31 (" divV) +

£ + div (pvV) = -g - 31 (, divV) +

d_

dx

du dv dy dx

f dv\ _ (3)

dy V dy/ '

Pcp + PcpVgradT = x^t + xVgradp + div (A gradT) + Qv, (4)

дт ^ дР

- + рсР^Т = х-

где х, у — декартовы координаты; V = (и, у) — скорость потока и ее проекции на оси х и у; р, р — плотность и давление; р — динамический коэффициент вязкости; ср — удельная теплоемкость при постоянном давлении; Т — температура; Л — коэффициент теплопроводности; Цу — мощность тепловыделения, обусловленного диссипативными процессами и нагревом от внешних источников энергии.

Уравнение сохранения энергии здесь записано в неконсервативной форме относительно температуры (в форме уравнения Фурье-Кирхгофа). Система уравнений (1)-(4) будет использоваться совместно с уравнением состояния идеального газа

р = РМт, (5)

где Д0 = 8,314 • 107эрг/моль-К — универсальная газовая постоянная; М — молекулярный вес газа, в общем случае МЕ = М^(р, Т). Функция х в правой части уравнения (4) представляет собой поправку на зависимость ср от давления.

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

де

Р^т + рVgrad е + р divV = Цу — divq, (6)

дЬ

где q = —Лgrad Т — вектор плотности потока теплоты за счет теплопроводности. В данной формулировке удельная внутренняя энергия е включает в себя все составляющие: энергию поступательного движения, вращательного, колебательного и электронного возбуждения атомных и молекулярных компонент газовой смеси. В уравнении (6) перейдем к энтальпии

к = е + р, (7)

р

для чего воспользуемся уравнением неразрывности (1). В результате чего получим

дк др Р + Р Vgrad к = — + Vgrad р + Цу — divq. (8)

дЬ дЬ

Если допустить, что имеется равновесие по всем внутренним степеням свободы частиц с энергией поступательного движения, то можно воспользоваться единой теплоемкостью, рассчитанной в приближении равновесной заселенности внутренних энергетических состояний, как, например, это сделано в работе [7]:

т

Н = 1 срйТ + Но, (9)

То

где Н0 — энтальпия образования вещества при Т = Т0.

Дальнейшее преобразование уравнения (8) будет основано на использовании свойств энтальпии как термодинамического потенциала:

Н = Н (р,Т), (10)

т. е.

дН\ /дН\ _ (дН\

dh= Ы тdp + Ы pdT = Ш T dp+c'iT- (11)

Это позволяет переписать уравнение (8) в виде

дТ iöh\ dp

+ pep VgradT + p^dhp) dp-

+ pV ( ^ ) gradp = + Vgradp + Qv - divq

föh\ \dP/'

или

дТ др

РСР-д£ + РсУ&а&Т = + Х^гаёр + - ё1УЯ, (12)

где

-р (дН)т- (13)

Воспользовавшись термодинамическими уравнениями Максвелла, получим

X =1 - Р

-■т- +1

дт \р)p p

= 1 -^^ . (14)

мД дТ )„

Отсюда следует, что в частном случае совершенного газа M^ = const, X =1.

Далее рассмотрим случай совершенного газа. Граничные условия задают невозмущенный набегающий поток

x = 0 : u = u0, v = 0, T = T0, p = p0, p = p0; (15) y = H : u = uo, v = 0, T = To, p = po, p = Po; ( )

и неотражающие граничные условия на выходе из расчетной области

^ du dv dT dp dx dx dx dx

(15б)

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

На поверхности обтекаемого тела задаются условия прилипания

y = 0: u = v = 0, T = Tw, дР = 0 .

дУ

(15в)

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

Рассмотрим преобразование системы уравнений (1)-(4) для ее интегрирования в произвольной криволинейной системе координат £ = £ (х, у), п = П (х,у). Уравнения (1)-(3) представим в следующем виде:

« Ж + = + ^ + Н; (16)

д£ дп

et +

f = 1

д£ дп p

pu pv

E = 1

F = 1

Ev =

1

J • Re

F =

-L V -

1

J Re

pU puU + CxP PvU + CvP

pV puV + nxP pvV + ny P

" 0

CxTxx + Cy Txy ixTxy + Cy Tyy

0

VxTxx + Vy Txy VxTxy + V'y Tyy

(17)

(18)

(19)

(20)

(21)

H = а

pv

- J

Txy puv

yJ •Re yJ

2р ' д (Чу v N + д f ПуЛ "

Re [д£ UjA

pv2

J

(22)

где U = £xu + £yv; V = Пх« + Пуv

Tra = J ^ ) 2

J

+

J / n

U)« V

J

v

+ а—;

yJ

Txy

ил /V *

") e + J n

Tyy = ^ 2

+ (M \

J )t+ V J Уп

U V v

j) « 4 J) n+V

и * = £у и + V * = Пу и + Пх^;

а = 0 — для плоского случая, а = 1 — для осесимметричного случая; £ = £ (х,у), п = П — криволинейные координаты, однозначно связанные с исходной декартовой системой координат; тхх, тху, туу — компоненты тензора вязких напря ений.

Уравнения (16)-(22) записаны в безразмерном виде, Яе = р0и0Ь0/р0 — число Рейнольдса, Ь0 — характерный размер.

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

J = det

£х £у Пх Пу

= £х % £y Vx

= £ = = д^ = дп

где£х = дх £у = д^ Пх = дх Пу = ду.

Напомним, что якобиан обратного преобразования 7-1 =

xx « XX п У« Уп

= x«Уп - У«Хп

позволяет рассчитать площадь в физическом пространстве для ячейки, заданной в расчетной области (для двумерного случая).

В преобразованиях систем координат весьма полезны соотношения

J=

1

пу £У

X« = -y, Хп = —y, У« =

J"« J' ~п J

Пх

" J:

Уп =

j

«

п

равнение сохранения энергии преобразуется к виду

Cp дТ U дТ V дТ

pJ Ж + pcp Jet + pcp J~ÖV

1 д_

Re Pr дС

X-

(tx + С2) дТ

дС

J

+

1 д Re Pr дС

X

(CxVx + Cy Vy ) дТ" дп

+ а

J

+

+

Re Pr yJ

1

Re Pr дп 1 д Re Pr дп X Г дТ

X

X

Ы + п

2) дТ

+

J дп

(CxVx + Cyпу ) дТ

J

дТ дп

Су дС + п^ ) + ~,

дС Qv J

+

(23)

где Рг = —-—0 — число Прандтля. Ао

Функция в правой части включает в себя следующие составляю -ие:

ЯУ _ Яга/^ Яр Яц (24)

Т = Т + Т + Т' ()

где Яг^ — объемная мощность тепловыделения, обусловленная радиационными процессами (в рассматриваемом случае полагается равной нулю); ЯР — объемная мощность тепловыделения, обусловленная работой сил давления; Яц — объемная мощность тепловыделения дисси-пативных процессов.

Два последних слагаемых представляются в виде

Qpl = (Uддр + _ui_.

J V JдC J дп J Cp0 То.

(25)

J = J 2

CxU

J

+

«

+

U Л (V*

uT) s H J.

п-xu

J

2

+2

+(w ^

J

i

J

+

U V v

j). 4 j) „ + vj

i

+ 2а —

+

u

Recp0 То

. (26)

Численный метод решения уравнений. Систему уравнений (16) для газодинамического этапа решения задачи представим в виде

f öEv дFy дt дС дп

д E„. + д F,

дС дп

+ H,

(27)

2

2

п

Tli li-i

M/2

ki+l/2

J

Pi

---Л i+1/2

Л i-1/2

Рис. 1. Расчетная ячейка

или

где

df dU f dE„ dV f dF„

--1---1--- +---1---

dt д£ d£ dn дп

д E„, + д F,

д£ дп

+ H, (28)

f = 1

p

pu pv

E- = J

0

£xP

£y P

F = — F J

0

ПхР Пу P

(29)

Для численного интегрирования (28) используются некоторые варианты конечно-разностных схем ЛИЗМ-типа [6]. Сначала рассмотрим исходный вариант ЛИЗМ-схемы.

Воспользуемся конечно-разностной сеткой

w

= {£,= 1, 2,..., N7; п, г = 1, 2,..., N1; т = ¿р+1 -

р = 0,1, 2,...}, (30)

и уравнение (28) проинтегрируем по площади заштрихованной ячейки (рис. 1).

Аппроксимация производных конвективных потоков через выде-ленну ячейку задается в следуем виде:

ди f = S«,j+1/2 - S«,j-1/2 ; д£ £j+1/2 - £j—1/2

ÖVf = ^Aj - S??

дП Пг+1/2 - Пг—1/2

i— 1/2,j

(31)

(32)

где

4«

Si,j+1/2

2 (MS+1/2 +MS+1/2I

pa/ J pua/J pva/J

+

1

+ 2 (MS+1/2 - |MUj+1/2|

pa/J pua/J pva/J

ij+1

S«,j-1/2

2 «-1/2 +|MU-1/2|

pa/ J pua/J pva/J

+

i,j-1

+ 2 (MU-1/2 - K-1/2I)

pa/J pua/J pva/J

, (34)

«J

sn

Si+1/2,j

= 1 (M

i+1/2j

+ IM,

V

i+1/2,j

pa/ J pua/J pva/J

+

ij

+ 1 (Mh/2,1 - |MV+1/2,j|)

pa/J pua/J pva/J

, (35)

i+1,j

Sn-1/2,j

2 (My-1/2,j + |MV-1/2,j|)

pa/ J pua/J pva/J

i-1,j

+ 1 (My-1/2,j - K-1Aj|)

+

pa/J pua/J pva/J

MS+1/2 = «+1/2) + + M+1/2) -,

MiU 1,

4 (mu+1)

при

«+1/2) + = 1 1

1 M + |MU|) при |Mg| > 1;

«+1/2)- =

- 4 (Miu+1 - 1)2

1 (Mj+1 - |MiUj+1|

при |Mg-+1|< 1, при |MiUj+1|> 1;

, (36)

(37)

(38)

(39)

MS-!/2 = «-1/2) + + -;

4 (Mg-1 + 1);

при

|MU_i| < 1,

«-1/2) + , 1

2 (M£._1 + при |Mj_1| > 1;

1

— (Mg - 1)2 при |M,

(MS_1/^ _ = { 1 4

1 (MU -|Mj|) при

|M j

I 4?

^ 1, > 1;

M U = —j.

ij a : a j

a ;j = 4 Y

,Pi ;j

Конечно-разностные формулы для §П+1/2 ^, ^, М- записыва-

ются аналогично.

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

(40)

(41)

(42)

(43)

(44)

(45)

(46)

д Ep SP;? Si;j + 1/2 SP;? Si;j_1/2

£¿ + 1/2 - &_1/2

д Fp Si+1/2;j SP,n Si_1/2;j

дп ni+1/2 - ni_1/2

р,?

S:

P+ =

0 0

P+ £x/J + P_ £*/J

. &/J . i;j . &/J -

i;j+1

Mj | < 1,

-Pij (1 + Mg) при |

при |Mj > 1;

Pi ;j

M- • + m. .

г} 1 г,} 1

Mu-

P_

^Pi ;j + 1 f1 - MS+1) пРи

при

mUm- IM.1

2

„ г,}+1 1 г,}+1| Pi ;j+1 M!j.+1

|Mj+1| < 1,

|Mj+1| > 1;

i,j-

0 0

P+ £*/J + P_ £*/J

. &/J . i_1;j . &/J -

i,j

1

-j (1 + mU-O

P+ =

2

2

Pi;j_1

1 Mj 1 + |Mj 1|

Mj 1

при при

|MU_1| < 1,

|Mj 1| > 1;

(47)

(48)

(49)

(50)

2рг,з (1 - Ми) при 1,

Рх Н 1 Ми.-Ши\ (51)

г^-^м^ при |м£|>1

В отдельных случаях (в зависимости от особенностей исследуемого течения) удается улучшить качество получаемых численных решений, применяя алгоритмы ограничения счетных потоков на гранях элементарных конечно-разностных объемов. Например, для ограничения

возможных численных осцилляций величин M-, My может быть ис-

JÜ Ъ/тУ

пользован MUSCL-алгоритм сглаживания [8].

Если ввести обозначение U для любой из величин M-, My, то

Ur = Uj+1 - 4 [(1 - к) Ai + (1 + к) A2]j+1 ; (52)

(A2)j = min mod [j - Uj); ß (U3 - U-i)]; (53)

(Ai)j = min mod [(Uj - Uj-i); ß (U+i - Uj)]; ( )

min mod (x, y) = sign (x) • max {0, min [x • sign (y); y • sign (x)]} ;

д 3 - к 1 (54)

1 < ß « w к =з-

По всей видимости, приоритет введения такой аппроксимации следует признать за работой Колгана [9], где она была сформулирована в виде

UL = Uj + - min mod (Uj - Uj-i; Uj+i - Uj); (55)

2

Ur = Uj+i - 2 min mod (Uj+i - Uj; U+ - Uj+i); (56)

a при (ab) > 0, |a| < |b|, min mod (a,b) = ^ b при (ab) > 0, |a| ^ |b|, (57)

0 при (ab) < 0.

Кроме возможного использования алгоритмов сглаживания, важным элементом вычислительной схемы является способ расчета скорости звука, по которой определяется число Маха в разложениях вида (37)-(51). Приведенные далее соотношения позволяют несколько сгладить решения на сильных разрывах:

a(i) = LP ,j+i/2 (58)

ai,j+i/2 = \ HP-, (58)

j ' V Pi >j+i/2

где

Pi,j+i/2 = 2 (Pi j + Pi,j+i) , Pi j+i/2 = 2 (Pi j + Pi,j+i ), (59)

или

,(2)

1

'ij+1/2 _ 2 (ai;j +

(60)

Конечно-разностная схема для уравнения сохранения энергии.

За основу возьмем уравнение сохранения энергии в следующем виде:

1 дТ - дТ V дТ

pcp j— + pcp + pcp = ,

(61)

где включает все слагаемые в правой части выражения (23).

Как уже отмечалось, в рассматриваемом случае уравнение сохранения энергии будет интегрироваться в неконсервативном виде. Имея в виду возмо ность использования излагаемого алгоритма для задач сильного радиационно-конвективного взаимодействия, для численного интегрирования (61) будет построена неявная конечно-разностная схема. Рассмотрим несколько альтернативных схем.

Схема Т-1. Получим конечно-разностное соотношение для расчетной ячейки, показанной на рис. 1. Ее объем равен

€.7 + 1/2 п. + 1/2

= ^ / ¿П = (£¿+1/2 - 6-1/2) (п^+1/2 - П — 1/2) = Рг&.

?j —1/2 Пj 1/2

(62)

Тогда, используя интегро-интерполяционный метод [3], получим

€. + 1/2 п. + 1/2

I? =

= P i

?j-1/2 nj —1/2 ( pcp— Т \

/pcp—дТ\ {—dn =

- ( pcp—т ) '

ij+1/2 V J /ij_1/2_

?j + 1/2

-Pi Ti;j / Щ (J ) d^

?j —1/2

(63)

или

= Pi

i,j

6+( ^r) + ^t) - 6+i pCpT

J

- b_( f T

i,j

i;j+1 PiTi;j

J

i,j_ 1

b+l ^ +

+ T

i;j + 1

- ^ J L 1 - « J

i,j

= pU Ti

м ij

r+i pcp \ l-l pcp \

M 7 - b- v 7 ■ v 7 ' ij-1 v 7 / i,j+i.

+

^ t __ j pCp \ r t ••+ f pCp

+ Ti,j+1bR —г - Ti,j- 1b + 7"

J

i,j+1

J

i,j 1

где

b± = 2 (Ui,j+i/2 ± |Ui,j+i/2|) , b± = 1 (Uij-1/2 ± |Ui,j-i/2|) , (65)

2

Ui,j±1/2 = 1 (UhJ +

(66)

Интегрирование в направлении ц производится аналогично. В итоге получается 5-точечная конечно-разностная схема

^+ Т^ + ^^ + Бг<3ТП+! - С^ТП+ + ^ = 0, (67)

где

Ai j —

a+ (pcp/7 )i -1,j a>

pi

+

Bi j =

a— (PCP/J)i+11, a+

P- Pi

i,j

pi

+ ; (68)

pi+pi

Ai,j =

b+ (PcP/J )i,j-i b—

qj

+

q3 qj

Bi,j =

b— (Pcp/j7+1 , b+

qj

+

+

qj+qj

(69)

= A^ + Bi,j + A^. В,7 Jp) 1;

V J J i,j т

(70)

= 1 (Vi+1/2,j ±

,j ± ' i+1/2,j

, a± = 1 (Vi_1/2,j ±|Vi_1/2,j|) ; (71)

Vi,j±1/2 = 2 (Vi,j + Vi±1,j)

(72)

Fi,j = T + Qdis,i,j + Qp,i,j + Qa,i,j + Qxy,i,j + Q Ji,j T

,i,j + Qp,i,j + Qa,i,j + Qxy,i,j + Qrad,i,j; (73)

5

Qdis,i,j

Ji,j Y"^ Q(k)

Re / v Qdis,i,j; k=1

(74)

Q(1) = 2

^x,i,j+1ui,j+1 £x,i,j— 1ui,j-1 ^ 1 +

Ji

i,j+1

+

Ji,j—1 / 2q,

nxJti+1,j Ui+1,j nx,i—1j ui—1,j\ 1

Ji+1,j Ji—1,j J 2pi

(75)

a

2

q(2) - 2

Ji

ij+1

+

J- 1 j 2qj

ny;i+1;j Vi+1;j ПУ , i _ 1 Vi __ 1 ;j

Ji

i+1 ; j

Ji_1 ; j

1

2Pi

Q(4)

^diS ; i ; j

Q.

(3) - 2a ( vij

dis ; i ,j " \ y T ; j Ji ; j

—i

i j+1

—i

i j_1

Ji

i j+1

Ji , j_1 / 2q

1

+

Vi

i+1 ,j

Ji

i+1 j

V %

i_1 ;j

Ji_1 ; j

1

2Pi

(77)

(78)

Q

(5)

diS ; i ;j

—i

i j+1

—i j 1 1

Ji

i j+1

Ji , j_1 / 2q

+

Vi

i+1,j

Ji

i+1,j

+

Vi_1,j Ji_1;j

1 + a Vi;j 2Pi yi;j Ji;j.

; (79)

Qp,i,j -

— Pi;j + 1 - Pi;j_1 + V Pi+1,j

2qj

i,j

2P

Pi_ j -L; (80)

Ji

i,j

ns+1

i+1 ;j+1

_Ts+1 + Ts+1 + T

Ti- 1 ;j 1 + 1;j+1 + T

s+1

-is+1

1;j+1 jLi+1;j_1

4Pi (

+ (ki,j + 1 - ki;j_1) (T^j - T_j) +

+

4Piij

(ki+1;j - ki_1;j ) (^Sj+1 - ^1)

4Pi^j

; (81)

Qa;i;j ^

Л.

i;j

1

Re Pr Ji,j

e (t s+1 _t s+M п (T s+1 _T s+1

sy,i,j VTi;j+1 Ti;j_1/ + 'У;i;j VTi+1,j Ti_1j

k = j - Re Pr

2qj

Лi;? (бхПх + Пу)

2Pi

i;j

Ji

i,j

11

- Ö (ah;i;j + ah;i±1;j ) ; ^ = X ( ^ ;i;j + ;

öft -

(п2 + n2)i,j

Re Pr Ji

I, Лi;j

bh - j

^ + ^ ij

Re Pr J,

i;j

(82)

(83)

(84)

(85)

P++ - п,+1 -п,; P_ - п,-П,_1; q+ - £¿+1; - 6-£ы; (86)

Qrad;i;j - (div qrad)

Lo

i;J

p0U0Cp ; 0 Ji ; j '

(87)

2

2

2

3

2

где — вектор плотности интегрального радиационного теплового потока.

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

а) Д+ = - , А- = - ¿7^,^-1;

б) если (Д+Д-) > 0, то Цу — не модифицируется;

в) если (Д+Д-) < 0, то = (1 - 2ш) + ш (¿/¿,^-1 + ¿¿¿,¿+1),

где ш < 0,5, ?7г,+ — значение сеточной функции, получаемой в результате численного интегрирования уравнений движения.

Схема Т-2. Это простейшая противотоковая схема 1-го порядка, при получении которой интегро-интерполяционным методом, функции рСрЦ/7 и рСр"/7 полагаются постоянными в пределах элементарного объема интегрирования, а их знак определяет способ расчета производных температуры. Пятиточечная конечно-разностная схема (67) в этом случае имеет следующие коэффициенты:

= — +—=—; Вг,+ = —+ + - +—; (88)

Р- Р- Рг Р++ Рг+Рг

= ^ + ^; = - ^ + 4+"; (89)

9? 9+ 9+9?

= + + + Aij + ; (90)

j т + Qdis ,i j

+ QP,i,i + Q

где

± __/ 1 I i\ _ pi,j cp,i,j

»c _ q (aij ± |aij; _

ij 1 / > T '

2

Остальные функции такие же, как в схеме Т-1. Схема Т-3. В этой схеме для расчета температурного поля используется ЛиБМ-алгоритм. Здесь так же, как и в схеме Т-1, уравнение сохранения энергии преобразуется к виду

дТ + д ( + д ( "Т)

РСРЖ + дё 1/ + дП 1 ) -

- Т| (рср 1)- Тдп (РсР 1) = ^. (92)

Интегрирование по объему расчетной ячейки дает

Tn+1 _ Tn

0п .сп. .Tj_3 +

т

1

+ —

Яз 1

+ — Рг

pCp UjTn+1) - (pcp Jtn+1

J J ij+1/2 V J / —1/2

m+1

pcpJT'+Л - (VpJT'

J / i+1/2,j V J J i—1/2,j

+

n+1

n+1 'i,j

11 U

Яз (л J

1 V V >

Рг и J)

ijj+1/2

i+1/2,j

U

PCP

J

рСр

J

i,j—1/2,

i—1/2jj

= . (93)

При расчете функций в круглых скобках используется АиБМ-аппроксимация. Например, для первого слагаемого

i U JjTп+1

яЛ J

rU

i,j+1/2

= ТП+1 (Mi,j+1/2

'г,3

+ ! ai,jpi,j ср,г,з

+

+ TT1, {MU

i,j+1 VMi,j+1/2 ±

qj J,j

— I ai,j + 1Pi,j + 1CP,i,j + 1

qj Ji

i,j+1

, (94)

MU+1/2

, а функции МУ^+1/2 рас-

где (<+1/2) = ^ (<+1/2

считываются по соотношению (37)

Аналогично представляются остальные слагаемые в выражении (93). Коэффициенты 5-точечной конечно-разностной схемы имеют следующий вид:

Ai,j =

+

Мг—1/2,3) / apcp \

Рг

v j л-1

j Р— Р

Вг .j =

,3 =

Вг ,j =

М V

М i+1/2,jj /apcp\ +

a

Рг

+

V J /г+1,3 Р+Рг

+

Mi,j—1/2) / apcp \ + b—

Яз ^ J 'ij—1 qj'

+

Му3+1/2) / apcp \ + b+

яз ^ j 7,з+1 qj+qj'

а остальные коэффициенты (Су, Fi¡j) и функции — такие же, как в схеме Т-1.

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

Результаты численного моделирования. Разработанный метод тестировался сравнением с результатами расчета параметров ударного слоя у сферы, обтекаемой сверхзвуковым потоком, представленными в работах [4,10]. В качестве примера такого сравнения на рис.2 показано распределение чисел аха у поверхности сфер радиусом Я = 66 см вязким теплопроводным совершенным газом с показателем адиабаты 7 =1,4 при следующих исходных данных: р« = 7,6 эрг/см3, рте = 2,82 • 10-8г/см3, МЕ = 29 г/моль, Тте = 297 К, V« = 1,035 х х 105см/с. Приведенные данные отвечают одной из наиболее грубых из использованных расчетных сеток (N1 = 55 х N7 = 101, где N1 — число узлов по нормальной к поверхности переменной (п)). Поэтому точность расчетов газодинамических функций у обтекаемой сферы, в особенности вблизи боковой и задней поверхности, является здесь весьма низкой. Тем не менее, приведенные на рис. 2 данные свидетель-

Рис. 2. Изолинии числа Маха при М =3, 7 = 1,4 вблизи поверхности сферы; кружками показаны результаты расчета [10] для положения ударной волны; х и у в см

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

На рис. 3 показано распределение скорости и (рис. 3, а), температуры (рис. 3, б) и давления (рис. 3, в) вдоль передней критической линии тока при различных числах Маха набегающего потока. Маркеры на линиях отвечают реальному местоположению значений сеточных функций. Для одного из вариантов (М = 10) сетка измельчалась в 2 раза. Последствия этого хорошо видны на рис. 3, а и рис. 3, б. При улучшении сетки наблюдается некоторое увеличение отхода ударной волны от поверхности тела и температуры в сжатом слое. На рис. 3, б показаны результаты расчета обтекания сферы при М = 20 и разных условиях теплообмена у поверхности сферы — охлаждаемой и теплоизолированной стенки. В случае теплоизолированной стенки толщина ударно-

- М=10

- М=10 (109x201)

- М=20

Рис. 3. Распределение продольной скорости и = —— (а), температуры Т, К (б) и

Р

давления Р = -— (в) вдоль передней критической линии тока от поверхно-

р™ У^

сти сферы (слева) до границы расчетной области вверх по потоку (справа) для разных чисел Маха; х в см

го слоя увеличивается примерно на 1 см. днако следует иметь в виду, что в целом температура сжатого слоя при Ы = 20 получается завышенной в силу использования предположения о совершенном газе. В действительности при таких скоростях движения уже заметно проявляются реальные свойства газа.

алее проанализирован результаты расчетов при тех е исходных данных, но с уменьшенным показателем адиабаты 7 = 1,2. В данной серии расчетов скорость набегающего потока V« изменялась в пределах 0,96 • 105 ... 3,2 • 105см/с, что обеспечивало диапазон чисел Маха Ы = 3 ... 10 и чисел Рейнольдса Яе = 975 ... 3250. Показатель адиабаты задавался отличным от обычно используемого значения 7 = 1,4 с целью получить температуру в ударном слое, не сильно отличающую -ся от той, которая получается для модели реального газа при больших скоростях. Таюке представляет интерес сравнение расчетных данных при двух указанных показателях адиабаты.

Расчеты проводились с использованием различных конечно-разностных сеток: N1 = 51 х N7 = 101; 101 х 201; 201 х 401 (рис.4). Из рисунка хорошо видна степень сгущения узлов сетки вблизи поверхности. ри построении этой сетки применялись аналитические методы, подробно изложенные в работе [11].

Сравнение результатов расчетов распределения давления и конвективных тепловых потоков по поверхности сферы от передней до задней критической линии тока показаны на рис. 5 и 6 (Ы = 3) и 10 соответственно. Из приведенных данных можно сделать заключение об удовлетворительной сходимости результатов при варьировании числа узлов сетки. Подтвержден также очевидный факт о большой чувствительности результатов расчета конвективных потоков от степени подробности сетки по нормали к поверхности. Расходимость результатов расчетов конвективных тепловых потоков, показанная на рис. 5, б, находится в прямом соответствии с качеством использованных сеток. Вблизи лобовой поверхности сферы конфигурации трех использован-нх сеток были близки, что дало практически одинаковые результат . В задней полусфере сетки заметно различались, что и явилось причиной расхождения плотностей конвективных тепловых потоков.

В расчетах использовались следующие соотношения для коэффициентов вязкости и теплопроводности:

-5т-г1,5 1 г .

р = 1,458 х 10-5T Л = р

110,4 + T' см • с' 7 Ro 1

(7 - 1) M Pro: где Pro = 0,72.

Рис. 4. Конечно-разностные сетки во всей расчетной области (а, в, д) и вблизи сферы (б, г, е):

а, б — сетка 51 х 101; в, г — 101 х 201; д, е — 201 х 401

Полное представление о поле течения у сферы при М = 3 можно получить, анализируя расчетные данные, показанные на рис. 7-11. Каждый из указанных рисунков состоит из двух частей. Первая часть (а) показывает распределение газодинамических функций во всей области течения, а вторая (б) — вблизи лобовой поверхности сферы. На этих рисунках показаны изолинии чисел Маха (рис. 7), давления (рис. 8), температуры (рис. 9), продольной и (рис. 10) и поперечной V (рис. 11) скоростей. Сравнивая указанные распределения, легко идентифицировать

вого потока, Вт/см2, (б) вдоль поверхности сферы для разных расчетных сеток при M = 3:

1 — сетка 51 х 101; 2 — 101 х 201; 3 — 201 х 401

основные элементы течения вблизи сферы: головную ударную (А) волну и сжатый слой (Б), звуковую линию у лобовой поверхности (В), веер волн разрежения (Г) и область донного течения (Д), свободный пограничный слой (Е), след (Ж), горло следа (З) и вязкое ядро следа (И).

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

р

120

110

100 90 80 70 60 50 40 30 20 10 0

ю1

q

10°

ю-1

10"2

10 0 50 100 150 5 СМ 200

б '

Рис. 6. Распределение давления, эрг/см3, (а) и плотности конвективного теплового потока, Вт/см2, (б) вдоль поверхности сферы для разных расчетных сеток при М = 10:

1 — сетка 51 х 101; 2 — 101 х 201; 3 — 201 х 401

ние газодинамических функций вдоль задней критической линии при М = 3 показаны на рис. 12. Из распределения числа Маха (рис. 12, а) видно, что в области возвратного движения М = 0,8. В точке с координатами х = 280 см М = 0. В этой точке продольная скорость меняет свое направление (см. рис. 12, б). Приведенные данные свидетельствуют о хорошей сходимости результатов распределений функций вдоль задней критической линии при измельчении сетки. На рис. 12, в и г

Рис. 9. Распределение температуры при М = 3 во всей расчетной области (а) и вблизи поверхности сферы (б); х и у в см

Рис. 10. Распределение осевой скорости (и/У(х) при М = 3 во всей расчетной области (а) и вблизи поверхности сферы (б); х и у в см

Рис. 11. Распределение радиальной скорости (у/У» ) при М = 3 во всей расчетной области (а) и вблизи поверхности сферы (б); х и у в см

показано распределение давления и температуры. Локальное разрежение сразу за сферой отвечает области вихревого движения. В окрестности задней критической линии тока температура возрастает до 356 К. Следствием локального сжатия газа на некотором расстоянии от сферы является возрастание температуры до 385 К. Указанное сжатие газа и его нагрев обусловлены сходящимся к оси симметрии течением в свободном пограничном слое.

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

Рис. 12. Распределение числа Маха (а), осевой скорости U = u/V^ (б), давления P = р/ (рж V2) (в) и температуры Т, K (г) вдоль задней критической линии тока от поверхности сферы (левая граница оси абсцисс) до конца расчетной области в следе (правая граница):

1 — сетка 51 х 101, 2 — 101 х 201, 3 — 201 х 401; х в см

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

Работа выполнена в рамках гранта РФФИ № 04-01-00237.

СПИСОК ЛИТЕРАТУРЫ

1. Surzhikov S. T. 2D CFD/RGD Model of Space Vehicles//Proc. of the Int. Workshop on Radiation of High Temperature Gases in Atmospheric Entry. October 2003, Lisbon, Portugal. European Space Agency, SP-533, 2003. - P. 95-102.

2. Su r z h i k o v S. T. Numerical Simulation of Heat Radiation Generated by Entering Space Vehicle // AIAA 2004-2379, 2004, 11 p.

3. Самарский А. А., Попов Ю. П. Разностные методы решения задач газовой динамики. - М.: Наука, 1980. - 352 с.

4. Белоцерковский О. М. Численное моделирование в механике сплошных сред.-М.: Наука, 1984.-519с.

5. М а р ч у к Г. И. Метода расщепления. - М.: Наука, 1988. - 263 с.

6. Edwards, J. R., Liou, M.-S. Low-Diffusion Flux-Splitting Methods for Flow at all Speeds // AIAA Journal. - 1998. - Vol. 36. - № 9. - P. 1610-1617.

7. Гурвич Л. В., В е й ц И. В., Медведев В. А. и др. Термодинамические свойства индивидуальных веществ. - М.: Наука, 1978. - 495 с.

8. I s s a, R. I., J a v a r e s h k i a n, M. H. Pressure-Based Compressible Calculation Method Utilizing Total Variation Diminishing Schemes // AIAA Journal. - 1998. -Vol. 36. - № 9. - P. 1652-1657.

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

10. БелоцерковскийО. М. Расчет осесимметричных тел с отошедшей ударной волной (расчетные формулы и таблицы полей течений). - М.: Вычислительный центр АН СССР, 1968.

11. Суржиков С. Т. Аналитические методы построения конечно-разностных сеток для расчета аэротермодинамики спускаемых космических аппаратов // Вестник МГТУ. Серия "Машиностроение". - 2004. - № 2. - С. 24-50.

Статья поступила в редакцию 29.03.2005 Сергей Тимофеевич Суржиков окончил МВТУ им. Н.Э. Баумана в 1975 г. Д-р физ.-мат наук, заведующий лабораторией "Радиационная газодинамика" Института проблем механики РАН, профессор кафедры "Теплофизика" МГТУ им. Н.Э. Баумана. Автор более 350 научных работ в области теплофизики и радиационной газодинамики.

S.T. Surzhikov graduated from the Bauman Moscow Higher Technical School in 1975. Dr. Sc(Phis.), Head of the Radiative Gas Dynamics Laboratory of the Institute for Problems in Mechanics Russian Academy of Sciences. Author of more than 350 publications in radiative gas dynamics and theory of heat and mass transfer.

В издательстве МГТУ им .Н.Э. Баумана в 2004 г. вышла в свет книга

Суржиков С.Т.

Тепловое излучение газов и плазмы. - М.: Изд-во МГТУ им. Н.Э.Баумана, 2004. - 544с.: 120 ил. (Компьютерные модели физической механики).

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

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

По вопросам приобретения обращаться по тел. 433-82-98; e-mail: [email protected]