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

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

УДК 533.6

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

АНАЛИТИЧЕСКИЕ МЕТОДЫ ПОСТРОЕНИЯ КОНЕЧНО-РАЗНОСТНЫХ СЕТОК ДЛЯ РАСЧЕТА АЭРОТЕРМОДИНАМИКИ СПУСКАЕМЫХ КОСМИЧЕСКИХ АППАРАТОВ

Приведены аналитические методы построения многоблочных конечно-разностных сеток, предназначенных для расчета обтекания и нагрева спускаемых космических аппаратов. Проанализированы преимущества и недостатки различных сеточных конфигураций. Представлены примеры расчетных сеток для трех форм космических аппаратов: Pathfinder (США), Mars Sample Return Orbiter (MSRO, Европейское космическое агентство) и сегментально-конического космического аппарата (Россия).

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

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

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

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

Построение конечно-разностной сетки вблизи КА сегментально-конической формы. Геометрия КА сегментально-конической фор -мы и двухблочная схема расчетной области приведены на рис. 1,а,б [1]. Первый блок расчетной области (Л^Б^С^ДО ограничен лобовой поверхностью КА (Л1Б1), внешней границей расчетной области со стороны набегающего невозмущенного потока (СхДх), осью симметрии (ДЛх) и конической поверхностью (БхСх), являющейся продолжением обратного конуса (Б2С2), образующего заднюю поверхность КА. Второй блок расчетной области ограничен осью симметрии 0Д2Е2, внешней границей, расположенной в невозмущенном газовом потоке (Л2^2), выходной границей газового потока (Е2^2) и границей Л2Б2С2Б2, состоящей из трех частей: Л2В2 — области газового потока, смежной к первому блоку расчетной сетки; Б2С2 и С2Д2 — ко -нической и сегментальной поверхностей КА. Направление изменения узлов расчетных сеток показано на рис. 1,а.

Для построения расчетной сетки заданы следующие исходные данные, характеризующие геометрию КА: Яп — радиус лобового сферического сегмента (радиус кривизны участка поверхности Л1Б1); Яь — радиус сферического сегмента донного затупления; 2< — угол раскрытия конической поверхности (< — угол между лучом Л2Б2С2 и осью симметрии); 2<2 — центральный угол лобовой поверхности (<2 — угол между лучом Х0Б2 и осью симметрии).

Заметим, что в практике экспериментальных и расчетно-теорети-ческих работ исследуется целое семейство сегментально-конических аппаратов, подобных тому, что показан на рис. 1. В настоящей работе проанализирована форма КА, задаваемая следующими значениями: Яп = 553, 8 см, Яь = 34 см, <1 = 35°, <2 = 22, 62°.

Заданы также параметры внешней границы расчетной области (см. рис. 1): Ь0 — осевая координата точки Е2; Я0 — радиальное расстояние точки

Положение КА внутри расчетной области задано только одной координатой Х0 — осевой координатой центра симметрии лобовой поверхности КА.

Используя заданные величины, можно рассчитать координаты точек Л1, Б1, С2, Д2, в, f:

= Х0 — ЯП ЯАх = 0;

ZBi = Xo - Rn cos RBl = Rn sin ;

i RBi D n

Zd = zbi + --; Rd = 0;

tg^i

. RB1 — RC2 t-, T-,

ZC2 = ZBi +--; Rc2 = Rb cos ^i;

tgP

Zf = ZC2 - Rb sin ^i; Rf = 0;

%в2 = %/ + Яь; Яв 2 = 0.

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

Рис. 2. Двухблочные расчетные сетки эллиптического (а), степенного (б) и параболического (в) типов

ординат или метода дискретных направлений (другое название этого метода — Ray-tracing method) [2]. Заметим также, что не последнее значение имеет эстетика расчетной сетки.

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

На рис. 2, а показана расчетная область, ограниченная эллиптической внешней границей. В трехмерном случае — это эллипсоид вращения. В этом случае уравнение внешней поверхности имеет вид

- + ^ + (Z - Lo)2 =1 (1) r>2 + -п2 + т2 1 (1)

R0 R0 L0

где ось z направлена вдоль оси симметрии КА; L0, R0 — полуоси эллипса в направлении z и r (см. рис. 1, б). Поскольку рассматривается

осесимметричная задача, то

2 2 I 2

г = х + у ,

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

r (z) = RH/ 1 - (z - Lo)2

L0

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

X - Хо = у - уо = г - го = £

Шх Шу

где х0, у0, го — координаты точки, из которой испускается луч; шХ,

—*

шу, <шг — направляющие косинусы единичного вектора О, характеризующего направление луча. Координаты указанной точки пересечения находятся при решении системы уравнений (1), (2), подробно исследуемой в аналитической геометрии [3]:

-В ± VВ2 - 4АС

¿1,2 =

2A

где

R

2

А = шХ + шу + Ш2

Д2

В = 2хоШж + 2уоШу + 2 (го - Ьо) 1;

Ьо

С = х2 + уо2 + (го - Ьо)2 ^ - Д*.

Из двух корней выбирается £ > 0:

Координаты точки пересечения с поверхностью определяются по уравнению (2)

х = хо + у = уо + Шу£, г = го +

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

Для построения сетки блока № 2 применен аналитический метод. На участке A2B2 расположение узлов сетки считается заданным, так как эта конфигурация уже определена для блока № 1. На участках поверхности B2C2 и C2D2 расположение точек задано в соответствии с функцией сжатия SBC и SCD (подробнее функции сжатия обсуждаются ниже). Кроме того, из общего числа точек NJ2, заданных для границ A2B2C2D2 и F2E2, число N11 определено для блока № 1. Оставшиеся (NJ2 — N11 + 1) точки распределяются на участках B2C2 и C2D2 в зависимости от геометрии задней поверхности КА. Таким образом, в рассматриваемом алгоритме построения расчетной сетки блока № 2 изначально считаются заданными координаты точек (rj, Zj) на границе A2 B2C2D2.

Из каждой точки (rj, Zj) до границы E2F2 проводится дуга эллипса, вычисляемая по следующим формулам:

r = j1 - > (3)

где

Ь; — . ± , Н; — — Ио.

'1 - М/Н) Г1

Соотношение (3) задает эллипс, проходящий через точки (г;, г;) и (Н;, Ь0) (см. рис. 1, б).

Заметим, что величину Н; можно выбирать по-другому, однако, если окажется, что Н; < г;, то указанный алгоритм приведет к ошибке.

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

Помимо расчетной сетки, ограниченной эллиптической поверхностью, представляет интерес класс расчетных сеток, ограниченных степенной зависимостью вида

г гт

— —-, (4)

Т г>т ' V '

Ь0 И0

где т > 2.

Как уже отмечалось ранее, уравнение любого луча, находящегося внутри расчетной области, определяется формулой (2), поэтому для определения координат точек пересечения луча с поверхностью необходимо решить систему уравнений (2) и (4), дополненную уравнением плоскости г — Ь0.

Решив уравнение (4) относительно г, получим нелинейную функцию следующего вида:

F (z) =

(

Ж0 + "<

z - Zo

"z

) + (yo + "y i-r)

Dm

- R-z = 0, (5)

Lo

для поиска корней которой можно использовать численный алгоритм. Простейший алгоритм представлен ниже.

При <шг > 0 имеем очевидное условие

г £ [го, Ьо] .

Если искомое решение г < Ьо, то оно отыскивается методом деления отрезка пополам (в качестве начального приближения выбирается [го, Ьо]). Проверим выполнение этого условия. Пусть г = Ьо, тогда из уравнения(2)получим

Ьо - го

t =

"z

что позволяет оценить неравенство

= \/(xo + "xt)2 + (yo + "yt)2 > Ro-

(6)

Если неравенство (6) справедливо, то луч пересекает поверхность (4). В противном случае луч пересекает плоскость г = Ьо. При < 0 имеем очевидное условие

г £ [0, го] .

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

Если = 0, то г = го, и расстояние от оси симметрии до поверхности _

п /г" г = до т . V Ьо

Поэтому для нахождения координат х и у надо решить квадратное уравнение

Г /12

(хо + Шж£)2 + (уо + Шу£)2 = До (го/Ьо)1/т ,

или

t2 + 2 ("xx + "yyo) t + X + y2 - R (zo/Lo)2/m = 0

2

2

r

так как в этом случае + — 1. Из двух корней уравнения следует выбрать положительный.

Заметим, что при ^ 0 изложенный алгоритм численного решения задачи о пересечении луча с внешней поверхностью может привести к большой погрешности расчета координат х и у. Поэтому в расчетах целесообразно применять алгоритм расчета для — 0 в некотором конечном диапазоне

Е [ — 8, +г] ,

где, например, г — 0, 01.

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

(г + шг Ь)т (г + £) ^

т

До Ьо

где шг — направляющий косинус луча по отношению к оси г.

Теперь рассмотрим построение расчетной сетки в 2-й области при использовании уравнения (4) для внешней границы. Пусть через точку (г;, г;) на поверхности КА проведена кривая, описываемая зависимостью (см. рис. 1, б):

г — гв гт Ьо — гв Щт

При построении сетки задаем изменение г^ по любому подходящему закону. Тогда

= rc ■ m , (7)

Lo - ZB

где

Tj

i-C — Rü"

TA

Rc — Rü—,

zb —

Zj - (ta/Rü)m L0

1 — (га/ДоГ '

Из рис. 1, б видно, что параметры Дс и гв задают координаты кривой (4), проходящей через точку (г;, г;).

Для расчетов по формуле (7) достаточно использовать однородную сетку г^ на отрезке г Е [г;, Ьо]. После того, как выполнен расчет координат (г;, г;) следует построить массив длин отрезков вдоль данной кривой:

li — Ii-i +

i + \/(т - Ti-i)2 + (zi - Zi-i)2; i — 2,..., N1; h — 0,

T

азатем, используя заданную функцию сжатия-растяжения, определить окончательные значения длин отрезков и соответствующие значения координат (т^- ). Конечно-разностная сетка у поверхности КА сегментально-конической формы при т = 3 показана на рис. 2, б.

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

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

Функции сжатия-растяжения узлов разностной сетки. Аналитические функции сжатия-растяжения позволяют создавать сгущение узлов разностной сетки в наперед заданных участках расчетной

Рис. 3. Двухблочная расчетная сетка степенного типа, используемая для аэротер модинамических расчетов

области. Такие функции создаются, как правило, для единичных отрезков, на которых вводятся переменные двух типов: по одной из них вводится однородная расчетная сетка, а по другой — неоднородная, с заданным аналитическим законом сгущения узлов сетки. Примеры различных функций сжатия-растяжения сетки можно найти в работах [4, 5].

Рассмотрим некоторые из алгоритмов. Пусть имеется некоторая криволинейная координата п, для которой необходимо создать неоднородную разностную сетку с заданным сгущением узлов (рис. 4). Пусть пА и пВ — начальная и конечная координаты, а ^ — Пз — ПА — длина соответствующего отрезка. Введем новую координату

Рис. 4. Криволинейная координатная линия в прямоугольной декартовой системе координат

П =

П - Па

Пв - Па

определенную на отрезке n е [0,1]. Введем однородную сетку

OJ

= {Vj, h = nj+i - nj = const; j =

От однородной сетки о к неоднородной сетке

о = {nj, h = nj+i - nj = const; j = 1, 2,...,NJ}

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

n = p j + (1 - P ) 1 -

{

th [Q (1 - jjj)] 1 th Q J

где Р и Q — числовые параметры, задающие степень неоднородности расчетной сетки.

В работе [5] проанализирована однопараметрическая формула вида

— 1п(1 + С^)

— 1п(1 + Q) ,

где Q — параметр, управляющий сгущением узлов.

Используя аналогичный подход, удается получить сгущение вблизи двух границ отрезка [5]:

Пз =

ln(1 + Pj) ln[l + Q (1 - j)]

И1,.........3

ln(1 + P)

ln (1 + Q)

}

Имеются также формулы для сгущения узлов сетки в окрестности промежуточной точки отрезка [6].

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

у = tg(в х),

где в — параметр сгущения сетки.

Рассмотрим построение сетки со сгущениями вблизи границ отрезка единичной длины [вА = 0, вв = 1]. Пусть

п п

хтт Т + и Хтах

2

2

где величины е1,2 = 0 так же, как и параметр сгущения сетки в, позволяют регулировать степень неоднородности сетки вблизи границ А и В (рис. 5). Если е1 = е2, то сгущение сетки вблизи обеих границ одинаково е, если е1 = е2 — разное.

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

{

шз = \ Уз = ymin + (j - 1) h; h =

Ут

Ут

NJ- 1

-; j = 1,2,...,NJ

}

Рис. 5. К построению неоднородной расчетной сетки с использованием тригонометрических функций

строится на отрезке (см. рис. 5)

ymin = tg (в x min У max = tg(e xmax)-

Тогда неоднородная сетка со сгущениями узлов вблизи границ A и B задается формулой

= «а + (sb - «а)—-mi^, Xj = — arctg(yj).

1

Xmax xmin ß

Если задать

п

xmin 0 и Xmax

тогда

ymin ymax tg(e Xmax)

и используя однородную сетку по х

~jj = ^ Xj _min , w nj _ i

ujj=j xj=xmin+(j -1) ^=xmaxj -x1min;j=i> 2,...,n j >

можно построить сгущение узлов в окрестности точки А:

(8)

в; — в а + ($в — в А) _УШ1" , У; — ^(в X;). (9)

утах ут1п

Если задать

п

хт1п ^ + 8 и хтах 0,

тогда

ут1п xm1n), утах 0,

используя однородную сетку (8) по х, можно построить сгущение узлов в окрестности точки В.

Еще одна модификация данного алгоритма позволяет вводить сгущение в окрестности промежуточной точки в* отрезка [вА, вв]. Величине в* поставим в соответствие величину у* — 0. Поэтому, если задать г — 0 и рассчитать

п

хтах — ^ — 8 Утах — Хтах^ то из выражения (9) найдем

в* вА / \

у* 0 ут1п + (утах ут1п)

вв — вА

и

__в* вА

ут1п утах ,

в - вв

после чего рассчитаем

_ 1

хтт в утт*

На отрезке [хт;п, хтах] задаем однородную сетку (8), после чего рассчитываем неоднородную сетку

Уз = tgв ,

а по выражению (9) — искомую шкалу —з.

Произвольные координаты п, используя сетки —з, преобразуют по формуле линейной интерполяции

Пз = ПА + (пв - пл) —-—.

вв - вА

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

Использование аналитических методов построения неоднородных сеток в двумерном пространстве. Простейшим и достаточно надежным методом, используемым в двумерном случае, является линейная интерполяция. Рассмотрим этот метод на примере построения сетки между осью симметрии Б2Е2 и внешней границей А2блока № 2 расчетной области (см. рис. 1, а). Пусть вдоль отрезков Б2Е2 и А2^2 заданы координаты (пбе)г и(п^)г, которые были рассчитаны с использованием любого из приведенных ранее методов построения неоднородной сетки. Тогда для каждой ]-й координатной линии

т,з = (ПАЕ)г + <Рз [(ПЯЕ)г - (плг,

где

ifj =

^,N.1 — 6,1

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

Х13, У1,з, (на границеА2 АО; ХNI,j, УN1,3, (на границе^^); Хг,1, Уг,1, (награницеА2^2>);

, Уг^1, (на границеДЕ2).

Рис. 6 (начало). Типы неоднородных сеток

Рис. 6 (окончание)

Тогда простейшая линейная интерполяция проводится по формулам:

xij = (1 - nij) xiNJ + Vijxi,i, (10)

yij = (1 - nij) Vi,NJ + VijVi,i, i = 1, 2, ...,NI; j = 1, 2,..., NJ. (11)

Сетка внутри расчетной области однозначно определена соотношениями (10) и (11).

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

Сетка «О»-типа для расчета обтекания КА Pathfinder. Геометрия КА Pathfinder (а также модель атмосферы Марса и параметры одной из траекторий входа) представлены в работах [7-9]. Для расчета обтекания КА Pathfinder применим расчетную сетку «О»-типа. Геометрия КА и расчетной области приведены на рис. 7. Расчетная область представляет собой полукруг диаметром TG = 2RD = 1000 см. Космический

Z, см

Рис. 7. Геометрия космического аппарата Pathfinder. Расчетная область "0"-типа

аппарат — это осесимметричное тело вращения (ABCDEF) с радиусом Миделя Rb = 130 см. Наветренная сторона КА состоит из конической поверхности BC, сопряженной со сферическим сегментом AB радиусом Rn = 66, 38 см. Угол полураствора конуса BC (по отношению к оси z) составляет = 70,19°. Подветренная сторона КА (DEF) представляет собой усеченную коническую поверхность DE. Угол по -лураствора данной поверхности = 46, 63°. Радиус задней круговой поверхности EF Ra = 28, 55 см. Конические поверхности BC и DE сопряжены круговым сегментом с радиусом Rc = 6, 62 см.

В данном случае (рис. 7) используется одноблочная расчетная сетка с пятью участками построения неоднородной сетки (см. рис. 7): ABST — область вблизи носового затупления; BCRS — область вблизи лобовой конической поверхности; CDPR — область вблизи сопряжения лобовой и задней конических поверхностей; DEHP — область за обратной конической поверхностью DE; EFGH — область следа.

Для расчета сетки на поверхности КА используются аналитические соотношения для каждого из пяти участков. Первый участок (AB):

Zi = zo - Rncosa; ri = Rnsina; a = (i - 1) hi;

i = 1, 2,...,M1; h1 =-1-(- -,

' ' ' 1 1 M1 - 1 \2 ^V '

где Ы1 — число узлов на первом участке, включая первый узел на оси симметрии (точка А) и последний узел (точка Б), являющийся одновременно первым узлом на втором участке; г0 — осевая координата центра кривизны носового сферического сегмента. Второй участок (ВС):

гг = гв + (г — 1) ^сов^ь г г = г в + (г — 1) г = 2, 3,...,М2;

, _ Ьво Т _ то — гв

п2 = Т7-7; ьво = —:-;

М2 — 1 эт<1

то = Яъ — Яо (1 — сов<1) ; тв = Япсов<ь

где М2 — число узлов на втором участке, включая точки Б и С. Заметим, что осевая координата точки С рассчитывается также по формуле

го = го — Япвт«р1 + Ьвосов<1.

Третий участок (СБ). Сначала рассчитываются координаты точки и, отвечающей миделеву сечению КА:

ги = го + Яссов<1 = го — Япвт<1 + £восов< + Яссов<1; ти = Яъ.

Тогда радиальная координата центра симметрии кругового сегмента СБ равна

ту = Яъ — Яс.

Координаты узлов сетки на поверхности этого сегмента рассчитываются по формулам:

гг = ги + Яс соваг; тг = ту + Яс втаг; аг = Н3 (г — 1) + <1, г = 2, 3,..., М3;

, <1 + <2

п3 =-,

3 Мз — 1 ,

где М3 — число узлов на третьем участке, включая граничные точки С и Б.

Координаты точки Б могут быть также рассчитаны по формулам:

гв = ги + Ясвт<2; тв = ту + Я^оя^.

Четвертый участок (БЕ). Осевая координата точки Е находится по формуле

Г в — те Г в — Яа

гЕ = гв +----= гв +----,

tg<2 tg<2

так как тЕ = Яа.

Координаты узлов сетки на конической поверхности

г^ = гв + (г - 1) Л,4^<2; Т = тв — (г — 1) , г = 2, 3,..., М4;

, ЬВЕ тп — те

п4 = Т7-Г; ьВЕ = —:-,

М4 — 1

где М4 — число узлов на четвертом участке, включая граничные точки Б и Е.

Пятый участок (ЕЕ). Координаты узлов сетки рассчитываются по формулам:

= гЕ; Ti = Те — (г — 1) Нь, г = 2, 3,...,МЬ; Нь = ГЕ ,

М5 — 1

где М5 — число узлов на пятом участке, включая граничные точки Е и Е.

Распределение точек на внешней границе ТО регулируется заданием долей полукруга ТО для отрезков ТБ, ТЯ, ТР, ТН, а также коэффициентами сжатия-растяжения сетки на каждом из отрезков.

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

Еще одним простым способом является использование одной или двух дополнительных поверхностей вблизи внешней границы области. Пример такой расчетной сетки показан на рис. 8. В этом случае координаты узлов двух вспомогательных поверхностей рассчитываются по формулам:

(2)

г( ) = + Ш2(гМ1}з — гМ1-2,1), 3 = 1, 2, ...,N3;

(2)

т( ) = Т N1-2,3 + Ш2 (т N1,3 — Т N1-2,3), 3 = 1, 2, ..., N3;

и

г^ = г* + ш^гшз — ¿32)), 3 = 1, 2, ..., N3; 3 = 3 + ^N1,3 — т32)), 3 = 1, 2,..., N3, где ш2 = 0, 5 и = 0, 8.

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

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

(2) (2) zNi-2j = zj ) ; rNI-2j = г) ) ;

ZNI-I,j = Z)(1); rNi-ij = j, j = 1, 2,NJ.

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

Трехблочная расчетная сетка «С»-типа для расчета КА MSRO. Форма космического аппарата MSRO (Mars Sample Return Orbiter [10, 11]) показана на рис. 10. Лобовая поверхность MSRO имеет три участка: AB — сферическое затупление радиусом Rn = 100 см (центр симметрии расположен на оси при z = z0), коническая поверхность BC с углом полураствора = 60°, скругление конической поверхности CD радиусом Rb = 15 см. Задняя часть КА MSRO представляет собой цилиндр радиусом R = 80 см и длиной 100 см. Представленная форма КА используется для расчетно-методических исследований обтекания космических аппаратов типа MSRO. Реальная форма проектируемого КА существенно более сложная, с несовпадающими осями симметрии аэродинамического щита и заднего цилиндра.

При построении расчетных сеток в данной работе используется метод аналитических отображений на трехблочной расчетной обла-

Рис. 9. Геометрия космического аппарата Pathfinder. Расчетная сетка "0"-типа 900

Zo

Z, см

Рис. 10. Геометрия космического аппарата MSRO и расчетная область "С"-типа

сти. Особенностью используемой системы сеток является изменение направления отсчета индексов между областями I и II и стыковка сеток в разных областях на смежных участках (БЕ) и (АВ)п, (БЕ)П и (АВ)ш. Внешняя граница расчетной области описывается единой степенной зависимостью вида (4).

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

Блок № 1 содержит три подобласти, построенные вблизи участков поверхности АВ, ВС и СБ. Расчетная сетка характеризуется следующими величинами:

= ,3 = 1, п ,% = 1,

% = 1,3 = 1 : точка А; г = 1,3 = NJl : точка Б; г = N11,3 = 1 : точка Е; г = N1-^,3 = NJ1 : точка Е} .

Первый участок (АВ):

г, = г0 — Rncosaj; г, = Я^та,; а, = (3 — 1) Н-, 3 = 2, 3,..., М-;

и ^

П- =

М- — 1

где Мх — число узлов на первом участке, включая точки А и В. Второй участок (ВС):

3 = гмг (3 — 1) К,2;

Гз = ГМ1 + (3 — 1) ПГ)2, 3 =2, 3,..., М2;

_ г1 — гь и _ гс — гь

пг,2 = "77 Т; пг,2 =

М2 — Г ',2 М2 — 1'

где М2 — число узлов на втором участке, включая точки В и С. Третий участок (СБ):

г, = гм2 + Яьсоя^х — Яь^ [<£- + Пз (3 — 1)] ;

г, = гм2 — Rьsin^l + Rьsin [^х + Пз (3 — 1)] , 3 = 2, 3,..., Мз;

= (Л/2)—* , 3 Мз — 1 '

где М3 — число узлов на третьем участке, включая точки С и Б.

Блок № 2 содержит две подобласти, построенные вблизи общей границы с блоком № 1 (АВ) и тыльной поверхностью аэродинамического щита (ВС). Расчетная сетка здесь характеризуется следующими величинами:

Шп = {£3,3 = 1, 2,...,^2; г/г, г = 1, ;

г = 1,3 = 1 : точка А; г = 1,3 = NJ2 : точка С;

г = N122,3 = 1 : точка Е; г = N122,3 = NJ2 : точка Б} .

Первый участок (АВ). Координаты точек на отрезке АВ совпадают с координатами узлов блока № 1 на участке БЕ. Здесь следует иметь в виду изменение порядка следования узлов в блоке № 2 по сравнению с блоком № 1.

Второй участок (ВС). Координаты узлов сетки, принадлежащих поверхности ВС, рассчитываются по формулам:

г3 = гА;

Тз = тв — (3 — 1) Ъ,2, 3 = NINI- + 1,..., N32;

, Тв — Тс

П,2. =

N3 — NI1

На цилиндрическом участке поверхности СБ строится сетка по индексу г. Используются формулы:

гг = гА + (г — 1)

тг = тс, г = 2, 3,...,ВД.

Блок № 3 содержит три подобласти, построенные вблизи общей границы с блоком № 2 (АВ) и донной поверхностью цилиндрического объема (ВС). На общей границе с блоком № 2 введены два участка, которые соответствуют двум участкам на левой границе второго блока. Сетка блока № 3 характеризуется следующими параметрами:

Шш = {£.3,3 = 1, 2,...,^з; Пг, г = 1, 2,...,NIз;

г = 1,3 = 1 : точка А; г = 1,3 = N1 : точка Е2;

г = 1,3 = N.122 : точка В; г = 1,3 = N3 : точка С;

г = N^,3 = N3, : точка Б; г = NI3,3 = NJ2 : точка Е3; г = NI3,3 = 1 : точка Е} .

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

Сетка, построенная с использованием представленных соотношений, функций сжатия и алгоритма генерации двумерных сеток, показана на рис. 11.

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

500

§ 400

с

300

200

100

0

0

250

500

750

1000

Z, см

Рис. 11. Геометрия космического аппарата MSRO и расчетная сетка "С"-типа

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

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

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

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

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

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

—*

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

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

Задача формулируется следующим образом. В локальной системе координат задан вектор

О = шх г I + Шу31 + ки

где шх, шу и <шг — направляющие косинусы вектора в локальной систе-

—* —* —*

ме координат с направляющими ортами г1, 31 и к1:

шх = cos<l ятв1; шу = ятв^ <шг = cosвl.

—*

Углы в1 и <1 отсчитываются от положительных направлений ортов г 1 и к.

Требуется найти направляющие косинусы этого вектора в лабораторной системе координат:

—* —>—>—>

О = шх г + шу 3 + шх к,

где шх, шу, шх — направляющие косинусы того же вектора в лабора-

—* —* —*

торной системе координат с направляющими ортами г, 3 и к.

Искомая связь устанавливается особенно просто в случае осевой

симметрии, когда единичные орты оси у лабораторной системы коор-

—* —*

динат~ и локальной системы координат 31 совпадают (рис. 12):

шх = ш%ятв + шхсояв; шу = шу; шх = ш%сояв — шхятв,

где в — угол между осями г и 2 или, что одно и то же, угол между локальной нормалью к поверхности и осью 2 лабораторной системы координат. Используя расчетную сетку, связанную с поверхностью КА, можно легко рассчитать косинус угла в:

_ Яа — Яь

^п I .

■у (Яа — Яь)2 + (2а — 2ь)2

i R i A Êf

. °ч

V, к / / / -

1 Л tz

► b l\

Z A Zb

Рис. 12. К пересчету угловых переменных в локальной и лабораторной системах координат

Формулы пересчета направляющих косинусов в трехмерном случае можно найти, например, в работах [12, 13].

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

рассчитываются по уравнению (2). Следующая задача как раз и состо-

—»

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

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

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

1. Румынский А. Н., Андреев В. М. идр. Укрупненный сценарий экспедиции для доставки марсианских образцов. Научно-технический отчет № 1549-1 по проекту МНТЦ № 1549. Центр исследований теплообмена ФГУП "ЦНИИМаш", 2000.

2. Filipskii M., Surzhikov S. Numerical Simulation of Radiation Heat Transfer in Plasma Generators. AIAA paper № 04-0988, Reno, NV, USA, 2004, 11 p.

3. Ильин В. А., Поздняк Э. Г. Аналитическая геометрия. - М.: Наука, 1981. - 232 с.

4. F l e t c h e r, C. A. J. Computational Techniques for Fluid Dynamics. 2. Specific Techniques for Different Flow Categories. Springer-Verlag, Berlin, Heidelberg, 1988.

5. Головачев Ю. П. Численное моделирование течений вязкого газа в ударном слое. - М.: Наука, 1996. - 376 с.

6. V i n o k u r, M. On one-dimensional stretching functions for finite-difference calculations. J. Comput. Physics, 1983, Vol. 50. No. 2. PP. 215-234.

7. M o s s J., et al. Journal of Spacecraft and Rockets, 1999, Vol. 36, No.3, pp.330-339.

8. C l a n c y R. T., Lee S. W., Gladstone G. R., Mc Millan W.W. and Rousch T. A New Model of Mars Atmospheric Dust Based on Analysis of Ultraviolet Though Infrared Observations for Mariner-9, Viking and Phobos. J. of Geophysical Research, 1995, Vol.100, No.E3, pp.5251-5263.

9. Milos F. S., Chen Y-K., Congdon W.M., Thornton J. M. Mars Pathfinder Entry Temperature Data, Aerothermal Heating, and Heatshield Material Response. J. of Spacecraft and Rockets, 1999, Vol.36, No.3, pp.380-391.

10. C h a r b o n n i e r, J. -M., F r a y s s e, H., V e r a n t, J. -L., T r a i n e a u, J. -C., Pot, Th., M a s s o n, A. Aerothermodynamics of the Mars PREMIER Orbiter in Aerocapture Configuration. Intern. Symp. on Atmospheric Reentry Vehicle and Systems, Acachon, France, March, 2001.

11. Verant, J. -L., Charbonnier, J. -M., B r o c, A., Diendonne, W., S p e l, M., Surzhikov, S. T., G r o m o v, V. G. Wake Flow Issues in Mars Aerocapture Phase. International Conference on the Methods of Aerophysical Research, Novosibirsk, 1-7 July, 2002, Vol. 3, pp.225-232.

12. Ольховский И. И. Курс теоретической механики для физиков. - М.: Изд-во МГУ, 1978, 574 с.

13. Ландау Л. Д., Лифшиц Е. М. Теоретическая физика. Механика. - М.: Наука, 1965.-203 с.

Статья поступила в редакцию 02.02.2004

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

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