Численный анализ высотной аэротермодинамики космических аппаратов тема автореферата и диссертации по механике, 01.02.05 ВАК РФ

На правах рукописи

Ващенков Павел Валерьевич

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

01.02.05 — Механика жидкости, газа и плазмы

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата технических наук

Новосибирск — 2012

005020135

Работа выполнена в Федеральном государственном бюджетном учреждении науки Институте теоретической и прикладной механики им. С.А. Хри-стиановича Сибирского отделения Российской академии наук

Научный руководитель: доктор физико-математических наук,

профессор

Иванов Михаил Самуилович

Официальные оппоненты: доктор технических наук,

Финченко Валерий Семенович

кандидат физико-математических наук Краус Евгений Иванович

Ведущая организация: Ракетно-космическая корпорация "Энергия"

имени С.П. Королёва

Защита состоится 20 апреля 2012 г. в 11:00 па заседании диссертационного совета Д.003.035.02 в Институте теоретической и прикладной механики им. С.А. Христиаиовича СО РАН по адресу: 630090, Новосибирск, ул. Институтская, 4/1.

С диссертацией можно ознакомиться в библиотеке ИТПМ СО РАН. Автореферат разослан " марта 2012 г.

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

Ученый секретарь диссертационного совета,

доктор технических наук

Засыпкин И.М.

Общая характеристика работы

Актуальность темы. Для успешного создания и использования космических аппаратов (КА) различного назначения необходимо детальное знание их аэротермодинамических характеристик (АДХ) вдоль всей траектории полета еще на этапе проектирования. Траектория К А проходит через зоны с различными режимами обтекания. Эти режимы характеризуются числом Кнудсена (Кп), определяемым как отношение средней длины свободного пробега молекул к характерному размеру КА. В орбитальном полете и на начальном участке траектории спуска с орбиты КА находится на большой высоте (более 150 км, Кп » 1), при свободномолекулярных условиях обтекания. Его аэротермодинамические характеристики определяются столкновениями молекул набегающего потока газа с поверхностью без учета межмолекулярных столкновений. Для изучения таких течений используются методы свободно-молекулярной газовой динамики. Обтекание КА на высотах ниже 60-70 км происходит в условиях сплошной среды. Числа Кнудсена в этих условиях достаточно малы (Кп < 1). Для исследования сплошносредных течений используются континуальные методы. Между этими предельными режимами КА проходит переходной режим обтекания, когда необходимо учитывать как столкновения молекул набегающего потока с поверхностью, так и межмолекулярные столкновения." Переходной режим характеризуется числами Кнудсена: 10_3 < Кп < 10. В этих условиях сплошносредные методы исследования неприменимы вследствие высокой разреженности и термохимической неравновесности газа. Экспериментальное моделирование высокоскоростных разреженных течений довольно проблематично, и в настоящее время практически единственным средством получения информации об аэродинамической обстановке около КА на больших высотах полета являются методы вычислительной аэродинамики. Для анализа аэротермодинамических характеристик КА в переходном режиме необходимо использовать кинетический подход (решение уравнения Больцмана или его моделей). Математическое моделирование аэротермодинамики КА в переходном режиме является весьма сложным и требует больших вычислительных и временных ресурсов. Использование таких трудоемких расчетных методов не всегда является целесообразным. Например, на начальном этапе проектирования К А необходимо исследовать его АДХ в широком диапазоне изменения праметров набегающего потока, углов атаки и в различных геометрических компоновках. Произвести эти расчетные исследования ресурсоемкими методами, основанными на решении уравнения Больцмана, за разумное время невозможно. Поэтому на начальном этапе проектирования КА для многопараметрических исследований АДХ применяются быстрые инженерные методы, а решения, полученные кинетическими методами, используются для оценки точности и проверки особых случаев. В связи с этим дальнейшее развитие и усовершенствование методов и средств

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

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

На защиту выносятся следующие разработки и результаты, составляющие научную новизну работы:

1. Алгоритмы балансировки загрузки процессоров для параллелизации метода прямого статистического моделирования Монте-Карло (метод ПСМ), основанные на минимизации межпроцессорного обмена;

2. Модификация локально-мостового метода для вычисления тепловых характеристик космических аппаратов в переходном режиме;

3. Результаты трехмерных расчетов аэротермодинамических характеристик КА "Клипер" в широком диапазоне чисел Кнудсена;

4. Результаты осесимметричных и трехмерных расчетов аэротермодинамических характеристик перспективной пилотируемой транспортной системы (ППТС) в гиперзвуковом потоке с учетом эффектов реального газа;

5. Результаты численного моделирования истечения струи из сопла двигателя ориентации в вакуум и формирования обратного течения;

6. Результаты численного моделирования трехмерного обтекания носовой части КА "Прогресс" после сброса створок головного обтекателя.

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

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

4

характеристики в пшерзвуковом режиме при многовариантных исследованиях модели КА. Этот метод включен в программный комплекс RuSat.

Программные комплексы SMILE и RuSat внедрены в Ракетно-космической корпорации "Энергия" и активно используются для исследования АДХ КА различного назначения (имеются акты о внедрении).

Исследован ряд прикладных задач аэродинамики разреженного газа и получена новая информация об:

- аэротермодинамических характеристиках КА "Клипер" и ППТС в широком диапазоне высот полета;

- особенностях течения около элементов КА, таких как разворот струи газа, истекающего из сопла двигателя управления в вакуум;

- распределении давления на конструкционных элементах и поверхности носовой части КА "Прогресс".

Эта информация необходима для проектирования и эксплуатации этих КА и содержится в их аэродинамических базах данных.

Апробация работы. Основные результаты диссертации докладывались и обсуждались на следующих конференциях и симпозиумах: 4-th European Symposium on Aerothermodynamics for Space Applications (Capua, Italy, 2001), 3-rd Atmospheric Reentry Vehicles & Systems (Arcachon, France, 2003), European Congress on Computational Methods in Applied Sciences and Engineering (Jyvaskyla, Finland, 2004), AIAA Thermophysics Conference (Portland, USA, 2004; San Francisco, USA, 2006; Honolulu, Hawaii, 2011), Всероссийский съезд по теоретической и прикладной механике ( Нижний Новгород, Россия, 2006), European Conference for Aerospace Sciences (Moscow, Russia, 2005; Brussels, Belgium, 2007), East-West High Speed Flowfields Conference (Beijing, China, 2005), International Conference on the Methods of Aerophysical Research (Novosibirsk, Russia, 2008), International Symposium on Rarefied Gas Dynamics (Monopoli, Italy, 2004; St.-Petersburg, Russia, 2006; Kyoto, Japan, 2008; Pacific Grove, USA, 2010)

Публикации. По теме диссертации опубликовано 20 печатных работ, список которых представлен в конце автореферата. 2 из них опубликованы в журналах из перечня ВАК. Кроме того, имеется свидетельство о государственной регистрации программы SMILE для ЭВМ.

Личный вклад автора. При выполнении работ по теме диссертации автором были разработаны модификации программных комплексов RuSat и

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

Структура и объем диссертации. Диссертация состоит из введения, четырех глав, заключения, списка цитируемой литературы из 48 источников и 64 рисунка. Полный объем диссертации 119 стр.

Содержание работы

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

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

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

Традиционно для параллелизации метода ПСМ используется декомпозиция расчетной области. При этом необходимо учитывать некоторые особенности метода: 1) вычислительная нагрузка на процессор зависит от количества частиц в ячейке, следовательно, количество частиц является индикатором загрузки; 2) в ходе расчета структура течения изменяется, в связи с этим изменяется количество частиц в различных зонах области; 3) значительная

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

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

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

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

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

Алгоритм с "делением области на равные части" является модификаци-

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

Работа алгоритма "жизнь" похожа на известную игру, моделирующую размножение колонии микроорганизмов. На начальном этапе в расчетной области выбирается по одной ячейке для каждого процессора. Эти ячейки будут являться центрами "жизни", и вокруг них будут наращиваться области ответственности процессоров по следующим шагам: 1) ищется область с минимальным количеством частиц; 2) в найденной области находится ячейка с максимальным потоком частиц через границу области; 3) ячейка, в которую направлен этот поток, присоединяется к области. Затем все шаги повторяются до тех пор, пока все ячейки не будут присоединены к зоне ответственности какого-либо процессора. Трудозатраты на применение этого алгоритма постоянны для заданного количества ячеек.

Оценка эффективности алгоритмов была проделана на примере поперечного обтекания цилиндра потоком газа с Кп=0.1 и М=5. Данное течение характеризуется значительной неоднородностью поля плотности в различных областях пространства (рис. 1).

Исследования проводились на вычислительном кластере МВС-15000ВМ Объединённого суперком-пыотерного центра (г. Москва) и на кластере МВС-1000 Сибирского суперкомпьютерного центра СО РАН (Новосибирск). Эффективность параллели-зации определялось по классической формуле: е = (Тх/М)/Тм, где Тх/М "идеальное" время расчета (время вычислений на одном процессоре поделенное на число процессоров), а Тм - реальное время вычислений на М процессорах. Полученная эффективность представлена в табл. 1.

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

ивврг

тж

Рис. 1. Поле плотности около цилиндра. Кп = 0.01, М = 5.

Таблица 1. Эффективность алгоритмов для разного числа процессоров.

Алгоритм Эффективность

211роц. I 411роц. | 81Іроц. | ШІроц. | 32ііроц. | Ь4ііроц.

Кластер МВС-15000

Вероятностный Жизнь Деление Подв. гран. 0.91 0.95 0.94 0.94 0.88 0.92 0.90 0.97 0.83 0.90 0.92 0.99 0.72 0.83 0.78 0.95 0.42 0.47 0.59 0.67 и.оУ 0.19 0.48 0.46

Кластер МЬС-Ши 0

Вероятностный Жизнь Деление Подв. гран. 0.95 1.03 0.98 1.03 0.93 0.98 0.95 0.99 0.81 0.97 0.83 0.83 0.80 0.98 0.86 0.85 0.04 0.80 0.74 0.83

ний, и предлагаемые алгоритмы способны их ускорить.

Вторая часть главы посвящена исследованиям АДХ различных моделей К А в широком диапазоне чисел Кнудсена локально-мостовым методом. Локалыю-мостовой метод позволяет быстро получить АДХ при проведении большого количества многовариаптных расчетов. Суть метода заключается в том, что для вычисления АДХ тела в переходном режиме используется интерполяция предельных значений, полученных в континуальном и свободно-молекулярном пределах. Значения АДХ вычисляются для каждой элементарной площадки, а затем интегрируются по всей поверхности обтекаемого тела:

дае.м.е,...) + с сот,4* ■ (1 - дае,м,е,...)),

с =

с^бБ, с = {сх,су,сг,тх,ту,гпг}.

Здесь 9 - угол между направлением потока и нормалью к элементарной площадке в данной точке; М - число Маха, Ие - число Рейнольдса, 5 - площадь поверхности КА. Функция ^ называется мостовой функцией. Существуют различные способы ее вычисления. В диссертации приводятся результаты, полученные с помощью функций, предложенных Коппенвалльнером, основанных на логарифмической зависимости мостовой функции от числа Кнудсена в невозмущенном потоке (в дальнейшем будем обозначать ее 1), а также с помощью полуэмпирической мостовой функции, основанной на числе Рейнольдса, предложенной Котовым, Лычкиным и др. (в дальнейшем 2).

В диссертации предлагается использовать локально-мостовой метод для вычисления коэффициента теплопередачи (<%) КА. Вычисление сд в континуальном режиме проводится на основе метода Ли и Нагаматсу. Коэффициент

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

Ch(s,e) = Chi

Здесь в - расстояние вдоль линии тока, 6 - угол между направлением потока и нормалью к элементарной площадке в данной точке, 7 - показатель адиабаты, с^о - коэффициент теплопередачи в точке торможения:

Чс = ±і,/5Еї і h^r»

v / 1V 7 i/Reoo^

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

Коэффициент теплопередачи на элементарной площадке в свободномоле-кулярном режиме вычисляется из классической теории свободномолекуляр-ных течений:

где

х (X) = е-*2 + spKX (1 + erf (х)),

Qe - коэффициент аккомодации энергии на стенке, 5оо - отношение скорости набегающего потока к наиболее вероятной скорости молекул, Tw - температура стенки, Too ~ температура набегающего потока, S^.e = SooCOs(Q). В работе предлагается использовать мостовую функцию вида:

2 V \ДКп2 6 VKn;

Если Кп0 < Кпт, используется мостовая функция і^д. В противном случае -

Значения Кпт = 0.3, ДКпі = 1.3 и ДКп2 = 1.4 были определены путем

сравнения с результатами моделирования методом ПСМ.

Рис. 2. Угловое распределение Си на сфере. рИс. 3. Коэффициент теплопередачи модели.

Эта методика использована в программном комплексе RuSat. Результаты расчетов, проведенных с его помощью, были сравнены с расчетами методом ПСМ. Так, например, зависимости коэффициента теплопередачи на поверхности сферы от угла при различных числах Рейнольдса, полученные в работе Dogra V.K., Wilmoth R.G., Moss J.N. "Aerothermodynamics of 1.6-m-Diameter Sphere in Hypersonic Rarefied Flow" (AIAA-91-0773) и вычисленные с помощью предложенной методики, показаны на рис. 2. Здесь угол откладывается от точки торможения. Как видно, тепловой поток и в континуальном, и в свободномолекулярном режимах определяется локальн-мостовым методом достаточно хорошо. Однако, в переходном режиме (Reoo = 3.8) локально-мостовой метод даёт завышенное примерно на 10% значение коэффициента теплопередачи на углах 20-60 градусов.

Кроме того, расчеты коэффициента теплопередачи локально-мостовым методом были проведены для цилиндрического тела с конической носовой частью. На рис. 3 представлены зависимости интегрального значения коэффициента теплопередачи этой модели с L = R (длина носика равна радиусу тела) от числа Рейнольдса в набегающем потоке. Здесь показаны кривые, полученные методом ПСМ (1), кривая, полученная из континуальной теории Ли и Нагаматсу (2), а также результаты локально-мостовой функции (3). Как видно, локально-мостовой метод дает гораздо лучшее приближение к результатам метода ПСМ при Re < 100, чем методика Ли и Нагаматсу.

Во второй главе приводятся результаты моделирования аэротермодинамики трех моделей КА "Клипер" (рис. 4) в широком диапазоне чисел Кнуд-сена. Две модели типа "Несущий корпус" отличались длиной аэродинамических щитков управления. Трехмерные расчеты на высотах от 80 до 150 км (7 • 10~4 < Кп < 100) были проведены методом ПСМ и локально-мостовым методом. На меньших высотах применение метода ПСМ требует недостижимых на сегодняшний день вычислительных затрат, и поэтому на высотах от

Рис. 4. Варианты геометрии КА "Клипер".

Рис. 5. Аэродинамические характеристики КА "Клипер" с короткими щитками при различных числах Кп. Угол атаки 0°. Кривые 1 и 2 - 1 и 2 локально-мостовые методы соответственно, кривая 3 - метод ПСМ.

50 до 80 км использовался только локально-мостовой метод.

На рис. 5 показаны зависимости аэродинамических коэффициентов сх, су, тг КА под нулевым углом атаки от числа Кнудсена в набегающем потоке. На рисунках нанесены результаты расчетов методом ПСМ, а также двумя локально-мостовыми методами, упомянутыми выше. Зависимости сх(Кп), полученные двумя локально-мостовыми методами, дают примерно одинаковое отличие от результатов метода ПСМ. Однако, зависимости нормальной силы ^(Кп) и коэффициента момента тангажа т2(Кп) различаются не только количественно, но и качественно. 2-й метод дает существенно лучшее совпадение с результатами метода ПСМ, чем 1-й. Он предсказывает даже такую особенность, как немонотонное поведение зависимости тг(Кп). Метод, основанный на логарифмической зависимости мостовой функции от числа Кнудсена неспособен показать такое поведение.

На рис. 6 показаны значения АДХ модели КА "Клипер" с удлиненными щитками под углом атаки 30° при различных углах отклонения балансировочных щитков. По этим кривым можно оценить эффективность щитков и возможность управления КА. Наиболее интересным здесь представляется график зависимости коэффициента момента тангажа от числа Кнудсена. Как видно, инженерные расчёты показали, что при числах Кнудсена Кп > 5 • 10"4 (высота более 80 км) коэффициент момента тангажа при любом отклонении балансировочного щитка отрицателен. Это означает, что КА невозможно будет сбалансировать на необходимом угле атаки 30° с помощью аэродина-

Рис. 6. Аэродинамические характеристики КА "Клипер" с удлиненными балансировочными щитками при различных числах Кнудсена. Угол атаки 30°. Кривые 1-3 - локально-мостовой метода для угла отклонения щитков 0, 10° и 20° соответственно. Точки 4-6 -метод ПСМ для угла отклонения щитков 0, 10° и 20° соответственно.

с„ м

мических органов управления. Лишь при полете на высотах менее 80 км с помощью балансировочного щитка можно получать как отрицательные, так и положительные моменты тангажа. Расчёты более точным методом ПСМ на высоте 80 км при различных углах отклонения балансировочного щитка показали, что погрешность локально-мостового метода приемлема, и полученные результаты корректны.

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

Третья глава посвящена исследованию аэротермодинамики перспективной пилотируемой транспортной системы (ППТС) (см. рис. 8) в диапазоне чисел Маха М=22-27 на высотах от 120 до 60 км. Исследовано влияние реальных свойств газа (возбуждение вращательной и колебательной энергии, диссоциации молекул) на АДХ КА. Также исследовано влияние эффектов реального газа на эффективность аэродинамических органов управления.

На высотах 75-110 км моделирование проводится методом ПСМ в осесим-метричной и в трехмерной постановке. Для задания химического состава набегающего потока была использована модель атмосферы ГОСТ-25645.154-90.

Рис. 7. Поле числа Маха и распределение коэффициента теплопередачи.

Рис. 8. Геометрия ППТС.

Рие. 9. Поле давления. Высота 75 км. Щиток отклонен на 30°. Химически реагирующий (вверху) и химически инертный газ (внизу).

На рис. 9 представлены поля давления около осесимметричной модели ППТС при отклоненном на 30° щитке на высоте 75 км. Вверху показан результат моделирования течения химически реагирующего газа, внизу - химически инертного. Отчетливо видно, что различия наблюдаются не только в расстоянии отхода головной ударной волны, но и во всей структуре течения около КА. В случае химически реагирующего газа температура потока за ударной волной значительно ниже. Распределение чисел Маха и конфигурация ударных воли в области течения также существенно зависит от учета химических реакций.

Интегральные значения ко-

Таблица 2. Коэффициент теплопередачи ППТС в

у эффициента теплопередачи пока-химически реагирующем и в химически инертном г,

газе|___ заны в таблице 2. Здесь можно

выделить резкое снижение влияния эффектов реального газа на тепловые нагрузки с увеличением высоты полета. На высоте 75 км значение, полученное без учета химических реакций, превышает значение в химически реагирующем газе в 3 раза. Как видно, влияние химических реакций проявляется на высотах менее 100 км. Это необходимо учитывать при проектировании теплозащиты КА. При этом коэффициент сопротивления зависит от наличия химических реакций очень слабо.

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

Для оценки эффективности аэродинамических органов управления бы-

Высота, км Реагир. газ Нереагир. газ

100 0.368 0.310

90 6.18-10"2 1.36 • 10" -1

85 3.38 • 10"2 9.47-10" -2

80 2.11 • Ю-2 6.22 ■ 10" -2

75 1.30-10~2 4.30 • 10" -2

ло проведено трехмерное моделирование течения около ППТС с щитками, отклоненными на 0, 20° и 30° на высотах от 80 до 110 км при углах атаки 0 и 40°. На рис. 10 показаны поле давления и поверхностное распределение коэффициента давления на высоте 80 км под углами атаки 0 и 40°. Из этих рисунков видно, что щиток управления вносит лишь слабое возмущение в поток при нулевом угле атаки, так как находится в следе за головной частью. При угле атаки 40° его влияние становится более существенным.

Рис. 10. Поле давления и поверхностное распределение коэффициента давления. Высота 80 км. Щиток отклонен на 30°.

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

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

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

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

Угол атаки а = 0

Угол атаки а = 40

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

Отработка этой методики производилась на экспериментальных данных , (Рот, 1971 г.). Характерные размеры сопла и течения: радиус горла сопла Д* = 2.55 мм, радиус выходного сечения сопла Де = 20.92 мм, температура торможения 7о = 300 К, давление торможения ро — 945 Па, число Рейнольд- 1 са Ие = 590. На рис. 11 показана исследуемая геометрия и поле плотности внутри и около сопла. Также показаны зоны течения, где использовались сплошносредные и кинетические методы. Внутри сопла течение моделировалось континуальным методом, прямоугольная зона на рисунке соответствует расчетной области, где использовался метод ПСМ. В качестве входных условий для метода ПСМ задавалась максвеловская функция распределения с параметрами, полученными из решения уравнений Навье-Стокса на части левой границы расчетной области для метода ПСМ, лежащей внутри сопла.

Для увеличения точности моделирования течения за пределами сопла была использована дополнительная подобласть в окрестности выходной кромки сопла. Подобласть для этого расчета представлена на рис. 12. Входные параметры на ее границе были взяты из предыдущего расчета методом ПСМ. Для того, чтобы оценить источники формирования возвратного потока, были проведены два расчета, отличающиеся тем, что в одном из них в моделируемую область был введен экран (рис. 12). Роль этого экрана заключалась в том, чтобы отсечь часть обратного потока, рождающегося в ближнем поле струи. Такой экран позво- ' ляет учесть вклад в обратный поток только молекул из пограничного слоя, разворачивающегося в окрестности кромки сопла.

На рис. 12 показаны поля плотности в области обратного течения, полученные без экрана и с установленным экраном около сопла. Как видно, изоли- | нии плотности над внешней стороной сопла в этих двух случаях существенно 1 различаются. Количественные данные представлены на рис. 13. Здесь показано радиальное распределение плотности за экраном, взятое на небольшом удалении от экрана. Введение экрана приводит к существенному уменьшению плотности. Это означает, что значительная часть обратного потока формируется не в пограничном слое около кромки сопла, а в периферийной части

Рис. 11. Поле плотности внутри сопла и в ближнем поле струи.

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

Во второй части этой главы проведено исследование полей давления около носовой части КА "Прогресс". Во время вывода КА на орбиту сам аппарат и его внешние конструктивные элементы закрыты обтекателем. При достижении определенной высоты, на которой скоростной напор становится достаточно слабым за счет уменьшения плотности, обтекатель сбрасывается. С экономической точки зрения желательно сбрасывать обтекатель как можно раньше. Однако при сбросе на недостаточно большой высоте аэродинамическая нагрузка на элементы КА велика, и существует риск разрушения конструкций. В результате расчетов были получены численные данные об аэродинамических нагрузках на носовую часть КА "Прогресс" (рис. 14). На рисунке отмечено положение двух датчиков полного давления (БЭШ и ОБШ), которыми оснащен реальный КА.

Исследовалось течение под углами атаки О и -8.5°. Параметры набегающего потока соответствовали полету на высотах 83 и 87 км (табл. 3) Трехмерное моделирование течения проводилось для двух моделей КА: без надстроек, и для полной геометрии с надстройками. На рис. 15 показано поле давления и линии тока около модели. В окрестности внешних элементов КА поле течения существенно возмущено. Из-за конструктивных особенностей датчики давления нельзя было разместить в невозмущенном потоке, поэтому для интерпретации результатов натурных измерений необходимо было определить параметры потока около их входных отверстий. В таблице 4 представлены полученные результаты трехмерных расчетов в окрестности датчиков.

Таблица 3. Параметры набегающего потока.

Н, км Мое т0о, к Рос, Па. Ыеоо КПоо

83 7.55 190 0.6 3950 2.3 • М-4*

87 7.69 190 0.3 2300 4.1 • 10~3

Рис. 12. Возвратное течение около кромки сопла: слева - без экрана, справа - со экраном.

Рис. 13. Радиальное распределение плотности за экраном.

Рис. 14. Модель полной геометрии носовой части КА "Прогресс" с надстройками.

41.7

45.8 50.0

54.2

58.3

62.5 66.7

70.6 75.0

79.2

83.3 87.5

91.7

95.8 100.0

Рис. 15. Поле давления и линии тока около геометрии носовой части КА "Прогресс" с надстройками.

Таблица 4. Параметры потока в окрестности датчиков давления.

Упрощенная конфигурация, а = О

Датчик Р/Роо Т/Тс» М ОЦ

DSN1 16.2 6.99 1.55 7.16

DSN2 14.1 6.76 1.94 7.16

Упрощенная конфигурация, а = 8.5°

DSN1 14.5 5.90 1.31 5.87

DSN2 11.5 6.61 1.93 3.27

Конфигурация с надстройками, а = 8.5°

DSN1 20.3 6.2 1.23 3.0

DSN2 19.6 7.7 1.59 2.96

Известно, что методика вычисления давления в набегающем потоке по соотношениям Рэлея в разреженном газе дает завышенное значение. Поэтому для нахождения полного давления, которое измеряет датчик в разреженном потоке необходимо использовать поправочный коэффициент. Для определения поправочного коэффициента в зависимости от степени разреженности потока были проведены дополнительные расчеты течения внутри датчиков DSN при указанных в таблице парметрах потока. Эта зависимость представлена на рис. 16. Здесь р = Rep = Индексы 1 и 2 обозначают величины перед и за

ударной волной, параметры которой посчитаны по формуле Рэнкина-Гюго-нио, d - внутренний диаметр трубки Пито.

С учетом этих поправок было посчитано полное давление, которое должны измерять датчики Пито, установленные на КА. Полученные результаты сравнивались с результатами, измеренными в полете (рис. 17). Здесь показано давление, измеренное датчиками, установленными на КА "Прогресс-М", запуск которого состоялся 24.01.2001. По оси абсцисс отложено время, которое К А находится в полете. Линиями отмечены результаты натурных измерений,

Рис. 16. Зависимость полного давле- Рис' 17' Сравнение результатов моделиро-п „ вания с результатами эксперимента. Запуск

ния от числа Реинольдса. опт

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

Основные выводы и результаты работы

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

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

- Исследовано влияние эффектов реального газа на аэротермодинамические характеристики перспективной пилотируемой транспортной системы (ППТС) в диапазоне высот от 60 до 150 км. Показано значительное (на высоте 75 км - в 3.5 раза) уменьшение тепловых потоков при учете химических реакций.

- Показано, что при истечении струи из сопла в вакуум формирование обратного течения происходит не только в результате разворота погранич-

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

- Получены распределения давления около носовой части КА "Прогресс" после сброса створок головного обтекателя. Сравнение с результатами натурного эксперимента показало достоверность результатов численного моделирования.

Список работ по теме диссертации

1. Кашковский А.В., Ващенков П.В., Иванов М.С. Программная система для расчета аэродинамики космических аппаратов // Теплофизика и аэромеханика. 2008. Т. 15. No. 1. С. 79-91.

2. Ващенков П.В., Кашковский А.В., Иванов М.С. Алгоритмы оптимизации вычислений методом ПСМ на параллельных вычислительных кластерах // Вычислительные методы и программирование. 2009. Т. 10. С. 290-299.

3. Свидетельство о государственной регистрации программы для

ЭВМ No 2008611184 "Программный комплекс SMILE "Statistical Modeling in Low-Density Environment" (статистическое моделирование в средах с низкой плотностью)" от 6 марта 2008 года. Правообладатель: ИТПМ СО РАН. Авторы: Иванов М.С., Маркелов Г.Н., Гимелынейн С.Ф., Кашковский А.В., Жукова Г.А., Бондарь Е.А., Ващенков П.В., Никифоров С.Б.

4. Krylov A.N., Kotov V.M., Tokarev V.A., Shcherbakov N.A., Khokhlov A.V., Ivanov M.S., Vashchenkov P.V., Kashkovsky A.V., Markelov G.N. Numerical Modelling and Experimental Data Analysis of the Flow near Spacecraft "Progress-M" Nose after the Head Fairing Jettisoning //4-th European Symposium on Aerothermodynamics for Space Vehicles, Italy, Capua, april, 2001, SP-487, 2002, P. 307-314.

5. Ivanov M.S., Vashchenkov P.V., Markelov G.N., Kashkovsky A.V., Krylov A.N. Numerical Simulation of High Altitude Aerodynamics of the "Progress" Spacecraft // Proc. of 3rd Atmospheric Reentry Vehicles & Systems Symposium, Arcachon, France, March 24-27, 2003, CD-ROM

6. Ivanov M.S., Vashchenkov P.V., Markelov G.N., Kashkovsky A.V., Krylov A.N. DSMC Study of the Near - Continuum Flow near the Nose Part of the Spacecraft "Progress-M": AIAA Paper 2004-2688, 2004.

7. Vashchenkov P., Kashkovsky A., Markelov G., Krylov A., Ivanov M. Numerical simulation of the high-altitude aerodynamics of the "Progres" spacecraft // Proc. of European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2004) / ed. Neittaanmaki P. etc., Jyvaskyla, Finland, 24-28 July, 2004 - V. 2. - Paper 542. - 15p.

8. Ivanov M.S., Kudryavtsev A.N., Markelov G.N., Vashchenkov P.V., Schmidt A.A. Numerical Study of Backflow for Nozzle Plumes Expanding into Vacuum : AIAA Paper 2004-2687, 2004.

9. Vashchenkov P., Kashkovsky A., Ivanov M.. Rarefied Aerodynamics of the Clipper Reentry Vehicle // 1-th European Conference for Aerospace Sciences (EUCASS), Moscow, Russia, 4-7 шоля 2005, CD-ROM

10. Vashchenkov P., Kashkovsky A., Ivanov M.. Numerical Analysis of Aerodynamics of Reentry Vehicles in Wide Range of Knudsen Numbers // Proc. of East-West High Speed Flowfields, Beijing, China, October 19-22, 2005, P.368-376

11. Vashchenkov P.V., Kudryavtsev A.N., Khotyanovsky D.V., Ivanov M.S. DSMC and Navier - Stokes Study of Backflow for Nozzle Plumes Expanding into Vacuum // 24-th International Symposium on Rarefied Gas Dynamics Monopoli, Bari, Italy / AIP Conference Proceedings 762, 2005, P. 355-360.

12. Vashchenkov P.V., Kashkovsky A.V., Ivanov M.S. Efficient Parallelization of the DSMC Method // 25-th International Symposium on Rarefied Gas Dynamics, Санкт-Петербург, 21-28 Июля 2006. / Rarefied Gas Dynamics. Proc. of 25-th Int. Sym. on RGD, Novosibirsk, Publ. House of the SB RAS, 2007, P. 533-538.

13. Vashchenkov P.V., Ivanov M.S., Krylov A.N, Khokhlov A.V. High-Altitude Aerodynamics of the Clipper Spacecraft // 25-th International Symposium on Rarefied Gas Dynamics, Санкт-Петербург, 21-28 Июля 2006. / Rarefied Gas Dynamics. Proc. of 25-th Int. Sym. on RGD, Novosibirsk, Publ. House of the SB RAS, 2007, P. 796-801.

14. Ivanov M., Vashchenkov P., Kashkovsky A. Numerical Modeling of High Altitude Aerodynamics of Reentry Vehicles // Proc. of 9th AIAA/ASME Joint Thermophysics and Heat Transfer Conference/ AIAA-2006-3801, San Francisco, California, USA, June 5-8, 2006. P. 1-12

15. Ivanov M.S., Kashkovsky A.V., Gimelshein S.F., Markelov G.N., Alexe-enko A.A., Bondar Ye.A., Zhukova G.A., Nikiforov S.B., Vashchenkov P.V.

SMILE System for 2D/3D DSMC computations // 25-th International Symposium on Rarefied Gas Dynamics, Санкт-Петербург, 21-28 Июля 2006. / Rarefied Gas Dynamics. Proc. of 25-th Int. Sym. on RGD, Novosibirsk, Publ. House of the SB RAS, 2007, P. 539-544.

16. Ivanov M., Vashchenkov P., Dyadkin A., Krylov A.. High-altitude Aerother-modynamics of the Clipper Spacecraft // Proc. of 2-nd European Conference for AeroSpace Sciences, Brussels, Belgium, 1-6 July 2007. P. 1-7, CD-ROM

17. Vashchenkov P., Kashkovsky A., Dyadkin A., Krylov A., Ivanov M.. Numerical Study of High-Altitude Aerodynamics of the Clipper Reentry Vehicle // Proc. of 46-th AIAA Aerospace Sciences Meeting and Exhibit / AIAA 2008-1188, 2008.

18. Vashchenkov P., Kashkovsky A., Ivanov M.. High-Altitude Aerodynamics of the Clipper Spacecraft // Proc. of 26-th International symposium on Rarefied Gas Dynamics, Kyoto, Japan, CD-ROM, ISBN:978-0-7354-0615-5.

19. Vashchenkov P., Ivanov M., Krylov A. Numerical Simulations of High-Altitude Aerothermodynamics of a Promising Spacecraft Model // Proc. of 27-th International Symposium on Rarefied Gas Dynamics, Pacific Grove, California, USA, 2010. P. 1337-1342.

20. Vashchenkov Pavel, Kashkovsky Alexandr, Krylov Andrey, Ivanov Mikhail Specific Features of Aerothermodynamics of a Promising Reentry Vehicle // Proc. of 42-nd AIAA Thermophysics Conference, Honolulu, Hawaii, 2011 / AIAA-2011-3128. P. 511-520.

21. Vashchenkov P.V., Kaskovsky A.V., Krylov A.N., Ivanov M.S.. Numerical Simulations of High-Altitude Aerothermodynamics of a Prospective Spacecraft Model // Proc. of the 7-th European Symposium on Aerothermodynamics, Brugge, Belgium, 2011: CD-ROM: SP-692 Aerothermodynamics, ISBN : 978-92-9221-256-8, 2011.

Ответственный за выпуск П.В. Ващенков Подписано в печать 15.03.2012

Формат бумаги 60x84/16, Усл. печ. л. 1.0, Уч.-изд. л. 1,0. Тираж 100 экз. Заказ №3

Отпечатано в ЗАО «ДокументСервис» 630090, Новосибирск-90, Институтская 4/1

61 12-5/2316

институт теоретической и прикладной механики

им. с.а. христиановича

На правах рукописи

Ващенков Павел Валерьевич

Численный анализ высотной аэротермодинамики космических аппаратов

01.02.05 - Механика жидкости, газа и плазмы

ДИССЕРТАЦИЯ

на соискание ученой степени кандидата технических наук

Научный руководитель

д. ф.-м. н., проф.

Иванов Михаил Самуилович

Новосибирск - 2012

Содержание

Введение ......................................................................4

Глава 1. Алгоритмы эффективной параллелизации метода ПСМ 23

1.1. Оценка эффективности параллелизации метода ПСМ......25

1.2. Эмулятор параллельных вычислений методом ПСМ ......30

1.3. Динамические алгоритмы......................32

1.4. Сравнение эффективности алгоритмов..............36

1.5. Локально-мостовой метод вычисления коэффициента теплопередачи ................................43

1.6. Оценка применимости локально-мостового метода........49

1.7. Выводы по главе...........................53

Глава 2. Моделирование обтекания космического аппарата "Клипер" .....................................55

2.1. Параметры набегающего потока и геометрическая модель ... 55

2.2. Параметры расчета методом ПСМ.................56

2.3. Результаты расчетов. Бескрылая модель..............57

2.4. Влияние эффектов реального газа на аэродинамические характеристики ...............................65

2.5. Результаты расчетов. Крылатая модель..............66

2.6. Выводы по главе...........................69

Глава 3. Моделирование обтекания перспективной пилотируемой транспортной системы......................71

3.1. Параметры набегающего потока и геометрическая модель ... 71

3.2. Результаты трехмерных расчетов методом ПСМ.........74

3.3. Учет отклонения щитков в осесимметричной постановке .... 78

3.4. Результаты осесимметричных расчетов..............80

3.5. Выводы по главе...........................88

Глава 4. Исследование особенностей течений разреженного газа около элементов космических аппаратов...........89

4.1. Истечение струи из сопла в вакуум................89

4.2. Течение около кромки сопла....................92

4.3. Моделирование обтекания носовой части КА "Прогресс" .... 94

4.4. Геометрия модели и параметры течения .............96

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

4.6. Выводы по главе...........................108

Заключение..................................111

Литература..................................112

Введение

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

Различные режимы обтекания характеризуются числом Кнудсена Кп = Х/Ь, где Л - средняя длина свободного пробега молекул, и Ь - характерный размер. Течение является континуальным, если число Кнудсена стремится к 0. При изучении таких течений можно пренебречь микроструктурой газа и использовать для расчета уравнения Эйлера или Навье-Стокса. При числах Кнудсена, стремящихся к бесконечности, режим течения можно рассматривать как свободномолекулярный. В этом случае столкновения молекул с поверхностью тела играют определяющую роль. При конечных числах Кнудсена необходимо также учитывать и столкновения молекул между собой. Такой режим течения называют переходным. Между континуальным и переходным режимом можно выделить "пограничную" область околоконтинуальных течений. Традиционно для расчета таких течений используются урав-

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

Определение границы применимости уравнений Навье-Стокса основано на анализе значений числа Кнудсена. Предполагается, что нарушение континуального описания происходит при числах Кнудсена более 10~2. Реально граница применимости континуального подхода зависит не только от значения числа Кнудсена, но и других факторов, например, геометрии обтекаемого тела. Для многих гиперзвуковых течений, представляющих практический интерес, присуща большая вариация параметров течения в окрестности обтекаемого тела. Это приводит к тому, что в некоторых областях применим континуальный подход, в то время как в других областях течения необходимо учитывать разреженность. Например, при гиперзвуковом обтекании затупленного тела в наветренной области газ сильно сжимается и нагревается, проходя через ударную волну. В этой области течение является континуальным, но в донной части и ближнем следе течение может стать разреженным. Здесь уравнения сплошной среды неприменимы и необходимо использовать уравнение Больцмана. Другим примером течения, где наблюдаются большие вариации параметров газа, является истечение газа в вакуум. Здесь течение является континуальным внутри сопла и около выходного сечения, переходным в ближнем поле струи и свободномолекулярным в дальнем поле.

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

где ф - параметр течения (например, плотность).

Другой критерий для определения границы применимости континуального описания был предложен Бердом для струйных течений

V 8 р йв '

где М - число Маха, р - плотность, 7 - отношение удельных теплоемкостей. В [1] показано, что континуальный подход становится неприменим при значении В > 0.05.

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

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

2. Метод прямого численного интегрирования уравнения Больцмана. Этот метод включает в себя два основных этапа - оценку интеграла столкновений с помощью метода Монте-Карло и интегрирование дифференциального уравнения. Основным недостатком этого подхода является существенная зависимость его трудоемкости от размерности задачи и, как следствие, весьма ограниченное использование для решения двухмерных и, тем более, трехмерных задач. Более существенным ограничением является то, что в настоящее

К ЩосаЛ = д

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

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

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

С уменьшением числа Кнудсена резко увеличивается время моделирования методом ПСМ. Стремление уменьшить время расчета для околоконти-

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

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

Попытки использования в методе ПСМ других сеток, разработанных для континуального подхода, не привели к положительному результату. Например, использование криволинейной (body-fitted) или монотонно-лагранжевой сеток (MLG) приводит к увеличению времени расчета, соответственно, в 2-i-lO раз и в 3-^4 раза. Анализ влияния сетки на структуру потока около затупленного тела показал, что наибольшее влияние оказывает размер ячейки по нормали к телу, а размер ячейки вдоль поверхности тела может быть порядка локальной длины пробега Л, что существенно уменьшает вычислительную стоимость моделирования, чем при использовании рекомендованных значений размера ячеек | [1]. Наиболее перспективным представляется использование прямоугольных многоуровневых сеток, которые позволяют обеспечить пространственное разрешение в зонах сильных градиентов параметров течения и сохранить высокую численную эффективность.

Ключевыми требованиями метода ПСМ, определяющими его ресурсо-емкость, являются следующие: 1) размер пространственной ячейки должен быть не больше локальной длины свободного пробега; 2) в каждой ячейке должно быть достаточно много модельных частиц. Эти два условия объясняют требовательность метода ПСМ к вычислительным ресурсам: при решении двумерных задач необходимое число ячеек и частиц для моделирования пропорционально квадрату плотности, для трехмерных задач число частиц

Модель Бесстолкновительное

дискретных Уравнение Больцмана уравнение

частиц Больцмана

Динамика сплошной среды

Уравнения Эйлера

Уравнения Навье-Стокса

Уравнения сохранения не образуют замкнутую систему

—I-1-1-1-1—

0.01 0.1 1 10 100 ->-00

Невязкий Локальное число Кнудсена Свободно-

предел молекулярный

предел

Рис. 1. Области применимости математических моделей описания течения газа (из [4])

меняется как куб плотности. Типичное количество модельных частиц для расчета обтекания КА с характерным размером порядка 1 м на высотах ниже 90 км методом ПСМ измеряется десятками и сотнями миллионов. Для расчета необходимо затратить порядка 105 процессоро-часов.

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

Расчет нестационарной стадии занимает довольно большую часть всего

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

Балансировка загрузки процессоров, применяемая в континуальных методах, основана на том, что время расчета каждого узла расчетной сетки примерно одинаково и не меняется со временем. Поэтому балансировка таких задач сводится к разбиению области на зоны с равным количеством узлов расчетной сетки, причем эти зоны не требуют перестроения. Такая техника разбиения также применяется и для метода ПСМ ([5, 6]). Однако, при использовании такого подхода в методе ПСМ невозможно получить равномерную загрузку процессоров на стадии установления вследствие существенного изменения числа частиц в каждой подобласти.

Специально разработанный для задач метода ПСМ алгоритм, в котором ячейки расчетной области распределяются между процессорами случайным образом [7], обеспечивает высокую равномерность загрузки процессоров, но совершенно не уделяет внимания межпроцессорному обмену.

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

Большинство используемых на практике параллельных алгоритмов балансировки загрузки процессоров для метода ПСМ [5, 9-14] обладают одним или несколькими из перечисленных недостатков:

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

- эффективны только для конкретной проблемы;

- применялись для маломасштабных задач;

- слишком сложны для реализации;

- Не учитывают время, затрачиваемое на передачу сообщений между процессорами;

- при увеличении числа процессоров увеличивается число сообщений.

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

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

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

провести анализ аэротермодинамики КА в широком диапазоне высот полета, параметров набегающего потока и в различных вариантах конфигурации КА.

Исторически первым пр