Построение метода Монте-Карло для решения задач высотной аэродинамики Текст научной статьи по специальности «Математика»

УДК 533.6.011.8

Зея Мьо Мьинт, А. Ю. Хлопков, Чжо Зин

Московский физико-технический институт (государственный университет)

Построение метода Монте-Карло для решения задач высотной аэродинамики

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

Ключевые слова: метод Монте-Карло, динамика разреженного газа, модели взаимодействия газа с поверхностью, модель Максвелла, модель Черчиньяни—Лампис— Лорда, космический аппарат, высотная аэродинамика.

1. Введение

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

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

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

Целью настоящей работы является исследование АДХ КА методом прямого статистического моделирования (Монте-Карло) в высокоскоростном потоке разреженного газа. В работе рассматриваются различные модели взаимодействия молекул газа с поверхностью и их влияние на АДХ.

2. Алгоритм метода прямого статистического моделирования (Монте-Карло)

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

/те •

Далее необходимо вычислить количество частиц, влетающих в область в единицу времени через все границы:

Щ = ^ I (£ ■ п3) (*, £)

)>0

где N3 — поток частиц через границу с номером п^ — единичный нормальный вектор. Вычисление сводится к известным интегралам от максвелловской функции, зависящим

т/ ь1/2 ь т

от скоростного отношения в = Утепте , пте = , ■

¿к! те

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

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

На третьем этапе по известным соотношениям, описанным выше, вычисляется скорость молекулы как случайная величина, распределенная в соответствии с функцией /те.

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

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

1 2

3

4

5

6

7

8 9

10

ввод данных;

определение номера части границы;

вычисление координат точки влета частицы в область;

вычисление скорости частицы;

вычисление координаты точки пересечения траектории частицы с поверхностью тела;

вычисление импульса и энергии, приносимых частицей;

вычисление скорости отраженной частицы;

вычисление импульса и энергии отраженной частицы;

выполнение пп. 4-7 до покидания молекулой расчетной области;

осреднение данных.

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

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

3. Модели взаимодействия молекул газа с поверхностью

Роль законов взаимодействия молекул с поверхностью проявляется тем сильнее, чем сильнее газ разрежен [1]. Граничными условиями для уравнения Больцмана являются условия, связывающие функцию распределения падающих и отраженных молекул. В модели Максвелла плотность распределения отраженных молекул имеет вид

fr(xw, ) = (1 - aT)fi (xw, - 2(£r ■ n)n) + aT nrexp (-hr, ■ n > 0,

и ядро рассеяния fl, 2] имеет вид

2h2

к(£i ^ с) = (1 - ^T)S & - 2 (£r ■ n) n] - Otexp [-hr£?] ■ & ■ n), hr = —.

Здесь полагается, что доля (1 - aT) молекул отражается зеркально, а остальная часть оу молекул — диффузно, параметр 0 < ат < 1 определяет коэффициент аккомодации касательной компоненты импульса aT = (PTi - PTr)/PTi.

Компоненты вектора скорости при диффузном отражении моделируются в локальной сферической системе координат, ось которой направлена вдоль вектора внешней нормали к поверхности, с помощью выражений [6]

|£r| = h-1/2\/- ln(aia2), cos в = д/од, ip = 2na4,

где ai, a.2-, од, — независимые случайные числа, равномерно распределенные в интервале (0, 1), в и (р — полярный и азимутальный углы.

Выражение для скорости отраженной молекулы с учетом неполной аккомодации по кинетической энергии имеет вид

|£r| = kh-1/2\/- ln(«i«2),

где к = \J(1 - (те )%ihr + (гЕ ■

Коэффициент аккомодации кинетической энергии определяется в виде

= Ei - Er = -Ei - Ew - h-1'

здесь Ew — энергия, которую уносили бы отраженные молекулы, если бы газ находился в равновесии со стенкой, т.е. когда Tr = Tw.

К. Черчиньяни и М. Лампис предложили феноменологическую модель (СЬ), которая также удовлетворяет принципу взаимности и является усовершенствованием максвеллов-ской модели [7]. Модель основана на введении двух параметров, которые представляют собой коэффициент аккомодации ап = оеп по кинетической энергии, связанной с нормальной компонентой скорости, и коэффициент аккомодации касательной компоненты импульса ат. Модель СЬ хорошо соответствует результатам лабораторных исследований с высокоскоростными молекулярными пучками [2]. Хотя сравнение ограничено лабораторными условиями, модель СЬ является теоретически обоснованной и относительно простой. Позднее появились модификации ядра рассеяния модели СЬ [8], однако они дают незначительное улучшение при сравнении с лабораторными экспериментами. В общем случае модель взаимодействия имеет несколько параметров произвольного физического смысла, которые позволяют добиться разумного согласия с результатами лабораторных исследований в некотором диапазоне условий. Универсальная модель должна использовать ядро рассеяния, полученное на основе физического эксперимента в широком диапазоне чисел Кнудсена и скоростей потока [2].

В модели СЬ ядро рассеяния для нормальной к поверхности компоненты скорости имеет

вид

к(£ni ^ inr) = 1о\2^1 - Ъ

(?п \

,гаг| sпг п—:-I exp

а„

О

Сп г + (1 - VnKn i

(7„

2ж

I0(x) = — I ехр(жcos -к J

здесь 1о — функция Бесселя первого рода, £п£пг — нормальная к поверхности компонента

т—1/2

скорости для падающей и отраженной молекул, отнесенная к пш .

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

К (tri ^ &г ) =

1

л/пат(2 - ат)

exp

(&r - (1 - От)Сп)2

(Гт (- - От)

касательная к поверхности компонента скорости для падающей и отражен-

здесь

молеку

Ядро рассеяния удовлетворяет принципу взаимности и условиям нормировки:

1—1/2

ной молекул, отнесенная к nw

l£J ¡м(£i)K(& ^ С) = l£nr| fM(C)KЫг ^ -£r), j К(& ^ )d£r = 1,

inr >0

здесь J'm — максвелловская плотность распределения.

Использованное преобразование расширяет СЬ-модель для учета обмена вращательной энергией между газом и поверхностью [8]. Модель в таком виде называется моделью Черчиньяни-Лампис-Лорда (CLL). Потом были предложены модификации модели [9] для учета обмена колебательной энергией и расширения диапазона состояний рассеянных молекул. Модель CLL в настоящее время получила широкое признание в работах [10-17].

o

4. Результаты расчетов и обсуждение

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

Геометрия тела представлена набором треугольников, к преимуществам такого представления относится простота описания формы тела и простота вычисления аэродинамических характеристик, основным недостатком является негладкость границы тела. Триангуляция поверхности тела и экспорт в формат STL (STereo Litography) были произведены с помощью комплекса автоматизированного проектирования SolidWorks на основании исходной модели, импортированной из формата IGS. Формат STL является достаточно простым и поддерживает твердотельность.

а)

б)

Рис. 1. Геометрия спускаемого аппарата (а) и крылатого космического аппарата (б)

Рис. 2. Зависимости Сх(а) для спускаемого аппарата (tw = 0.04)

Значения параметров: температурный фактор = /Т^ = 0.04, 0.1; скоростное отношение 5 = 20; коэффициенты аккомодации тангенциального импульса и нормальной энергии ат, ап = 0.5, 0.75, 1. Расчет проводился с использованием 5 • 106 частиц.

На рис. 2-4 представлены зависимости коэффициентов силы сопротивления Сх, подъемной силы Су, момента тангажа шг от угла атаки а от 0° до 90° для спускаемого аппарата (рис. 1а). При уменьшении коэффициента аккомодации ат от 1 до 0.5 величина Сх увеличивается при 0° < а < 55° и значение достигает до 2.72 при а = 0°. Коэффициент Су снижается в несколько раз по модулю при уменьшении ат от 1 до 0.5 при 0° < а < 80°. Зависимость шг(а) объясняет то, что при ионижении ат чувствительно увеличивается на всех диапазонов угла атаки.

Рис. 3. Зависимости Су(а) для спускаемого аппарата (tw = 0.04)

тг

ч\ Г

■As—-JS—ät. Метод Монте-Карло (Maxwell, <тт = 0.5) ^—+—+ Ие-сд Монте-Карло (CLL, с, = 0.5) Ф——❖ Метод Монте-Карло (с, = 1) -»-X Метод Ньютона OL

1 1 1

10 30 50 70 00

О 20 40 ВО ВО 1ВВ

Рис. 4. Зависимости mz(а) для спускаемого аппарата (tw = 0.04)

На рис. 5-7 представлены зависимости коэффициентов силы сопротивления Сх, подъемной силы Су, момента тангажа mz от угла атаки а от -90° до +90° для крылатого космического аппарата (рис. 16). При уменьшении оу от 1 до 0.5 вели чина Сх снижается до 1.85 при -55° < а < 55°, и при уменьшении оу от 1 до 0.75 величина Сх снижается до 1.74 при -55° < а < 55°. В рамках модели Максвелла при больших по модулю углах атаки зеркально отраженные молекулы повышают величину Сх, чего не наблюдается в рамках модели CLL. При уменьшении оу от 1 до 0.5 величина Сх увеличивается до 2.64 при а = +90°. Коэффициент Су снижается в несколько раз по модулю при уменьшении ат от 1 до 0.5, 0.75. Зависимость mz(а) объясняет то, что при понижении ат чувствительно

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

Рис. 5. Зависимости Сх(а) для крылатого космического аппарата =0.1)

Рис. 6. Зависимости Су (а) для крылатого космического аппарата =0.1)

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

При (Гт, ип = 0, 1 Сх и Су значительно совпадали. Но при малых углах атаки отра-

жснныс молекулы слабо отклоняются от первоначального направления и поэтому вносят малый вклад в сопротивление тонкого тела. При дальнейшем увеличении vivía атаки ситуация изменяется: зеркально отраженные молекулы передают поверхности тонкого тела больший импульс, чем диффузно рассеянные от стенки молекулы. Поэтому максимум сопротивления при больших углах атаки соответствует случаю ат = Ос полной аккомодацией ит = 1. Отметим, что близость результатов, полученных с помощью моделей Максвелла и CLL, отмечалась ранее в работе [13] для тел с высокими коэффициентами аккомодации поверхности, что позволяло достигнуть лучшего согласования с результатами эксперимента в аэродинамической трубе [17, 18].

-0-О Метод Монте-Карло (Maxwell, 'Гт - Li :

it.— ^-* Метод Монте-Карло (Maxwell, cx-ij.75J

□-В-□ Метод Монте-Карло (Maxwell - 1)

+-+--s- Метод Монте-Карло (CLL, :т- 0 :7

О-©-Q Метод Монте-Карло (ULL ах - 0,75)

ir-if-i* Метод Ньютона

-30 -75 -Б0 -45 -30 йЩ 0 15 30 45 Б0 75 90

Рис. 7. Зависимости mz(а) для крылатого космического аппарата (tw = 0.1)

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

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

Литература

1. Коган, М.Н. Динамика разреженного газа. М. : Наука, 1967.

2. Бараицев Р.Г. Взаимодействие разреженных газов с обтекаемыми поверхностями. М. : Наука, 1975.

3. Хлопков Ю.И. Статистическое моделирование в вычислительной аэродинамике. М. : МФТИ, 2006.

4. Белоцерковский О.М., Хлопков Ю.И. Методы Монте-Карло в механике жидкости и газа. М. : ООО «Азбука 2000», 2008.

5. Хлопков Ю.И., Звя М.М., Хлопков А.Ю., Т1жо 3. Методы Монте-Карло для определения аэротермодинами ческих характеристик гиперзвуковых воздушно-косми ческих

систем // Materials digest of LI International Research and Practice Conference «Physical, Mathematical and Chemical Sciences: Theoretical Trends and Applied Studies». — London: IASHE. - 2013. - P. 41-44.

6. Bird G.A. Molecular Gas Dynamics and the Direct Simulation of Gas Flows. — Oxford: Clarendon Press, 1994.

7. Cercignani C., Lam,pis M. Kinetic Models for Gas-Surface Interactions // Transport Theory and Statistical Physics. - 1971. - V. 1, N 2. - P. 101-114.

8. Lord R.G. Application of the Cercignani-Lampis Scattering Kernel to Direct Simulation placeMonte Carlo Calculations // Proc. of 17th Int. Svmp. on Rarefied Gas Dynamics. — 1991. - P. 1427-1433.

9. Lord R.G. Some Further Extensions of the Cercignani-Lampis Gas-Surface Interaction Model 11 Phvs. Fluids. - 1995. - V. 7, N 5. - P. 1159-1161.

10. Ketsdever A.D., Muntz. E.P. Gas-Surface Interaction Model Influence on Predicted Performance of Microelectromechanical System Resistojet // Journal of Thermophvsics and Heat Transfer. - 2001. - V. 15, N 3. - P. 302-307.

11. Utah S. and Arai H. Monte Carlo Simulation of Reentry Flows Based Upon a Three-Temperature Model // Proc. of 23rd Int. Svmp. on Space Technology and Science. — 2002. _ у. i. _ p. 1209-1214.

12. Santos W.F.N. Gas-Surface Interaction Effect on Round Leading Edge Aerothermodvnamics 11 Brazilian Journal of Physics. - 2007. - V. 37, N 2 A.

13. Padilla J.F. Assessment of Gas-Surface Interaction Models for Computation of Rarefied Hypersonic Flows // Ph.D. Dissertation. — University of Michigan, 2008.

14. Wadsworth StateD.C., Van Glider D.B., Dogra V.K. Gas-Surface Interaction Model Evaluation for DSMC Applications // Proc. of 23rd Int. Svmp. on Rarefied Gas Dynamics.

- 2003. - P. 965-972.

15. Воронин И.В., Мьинт З.М. Влияние особенностей взаимодействия газа с поверхностью на аэродинамические характеристики космического аппарата // Вестник МАИ. — 2010.

- Т. 17, № 3. - С. 59-67.

16. Зея М.М., Хлопков А.Ю., Чжо 3. Основные подходы к построению методов Монте-Карло в вычислительной аэродинамике // Труды МАИ. — 2011. № 42. — 17 с.

17. Padilla J.F., Boyd I.D. Assessment of Gas-Surface Interaction Models in DSMC Analysis of Rarefied Hypersonic Flow // AIAA Paper 2007-3891. - 2007.

18. Kussoy M.I., Stewart D.A., and Horstman C.C. Hypersonic Rarefied Flow over Sharp Slender Cones 11 NASA Technical Note. 1972. - D-6689.

Поступим в редакцию 12.10.2013.