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

УЧЕНЫЕ ЗАПИСКИ ЦАГИ То м IX 197 8

№ 4

УДК 533.6.011.8

ОБТЕКАНИЕ ОСЕСИММЕТРИЧНЫХ ТЕЛ ГИПЕРЗВУКОВЫМ ПОТОКОМ РАЗРЕЖЕННОГО ГАЗА

Ю. И. Хлопков

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

Интерес к исследованию течений для переходного режима в окрестности осесимметричных тел объясняется двумя причинами: с одной стороны, простая форма обтекаемого тела упрощает проведение расчетных исследований, с другой стороны, удается реализовать эффекты трехмерного обтекания реального летательного аппарата в простейшей форме. Этой проблеме посвящено большое количество работ, однако в основном они касаются двух крайних режимов обтекания: свободно молекулярного и режима сплошной среды. Для вычисления же аэродинамических характеристик на переходном режиме обтекания, как правило (особенно это касается конуса), пользуются приближенными полуэмпирическими методами. Исключение составляют методы, изложенные в работах {1, 2]. В работе [2] методом прямого статистического моделирования проводится расчет практически всех тел простой формы. Однако приведенные в этой работе данные не позволяют достаточно полно судить об аэродинамической картине обтекания осесимметричных тел. Трудность решения подобного сорта задач заключается в сложности решения уравнения Больцмана. Кроме обычных трудностей, возникающих при решении кинетического уравнения, при обтекании осесимметричных тел появляется дополнительная трудность, связанная с размерностью задачи. В данной работе методом Монте-Карло решается модельное кинетическое уравнение. Используемый метод [3] уже применялся для решения осесимметричных задач, но в данной работе проведено его усовершенствование.

1. Общие свойства метода. Рассмотрим кинетическое уравнение

<*//<# = /(/). (1) Будем решать его методом последовательных приближений

/о= /о>

и =/« + /(/о)Д^;

(2)

А*.

Тогда, считая, что схема установления (2) сходится, решение кинетического уравнения (1) можно представить в виде

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

Ф = (/> Ф)»

где ф — некоторые микроскопические характеристики газа.

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

™ N1 Р1(Ы ’ к = 1

где р — плотность вероятности, с которой распределены частицы;

I — номер ячейки в физическом пространстве.

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

«£<*• <4>

где Дл: — ячейка физического пространства.

Тогда для оценки погрешности при вычислении макропараметров за п шагов можно воспользоваться известной формулой [4]

8Фл < 8Фо + я д* О (Д/, Дх). (5)

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

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

/ (/) = .(/<)_/), (6)

где V — частота столкновений,

/° — равновесная функция распределения, то кинетическое уравнение (1) допускает запись в моментной форме

Таким образом, из схемы последовательных приближений (2), записанной для модельного уравнения, следует простая процедура счета

fn+l = /л (1 ~ vn bt) + /° V„ Д*.

При увеличении шага по времени на Дt пробная частица с вероятностью, пропорциональной (1—не меняет функции распределения

/(1 + Дг, х + Ш, ?)=/(*, a:, S)[l — v(jc) Д^];

с вероятностью, пропорциональной v„At, меняет функцию распределения на равновесную:

х + Ш, х, %)ч{х)Ы.

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

хп+1 =•*;„ + №.

С вероятностью 1—чМ частица не изменяет скорости, а с вероятностью vA/ ее скорость выбирается в соответствии с равновесной функцией распределения. Чтобы удовлетворить условию устойчивости (4), шаг по времени выбирается автоматически. Поле течения разбивается на ячейки, и шаг по времени выбирается из условия, при котором траектория частицы не пересекает границ ячейки. Накопление и вычисление макропараметров осуществляются обычным статистическим способом (3).

Интересно отметить, что при таком подходе к обоснованию метода [3] становится ясно, что описанная методика может быть использована не только для модельного уравнения Крука (6), но и для других форм модельных уравнений, например, для модели Халвея [6] и для аппроксимирующих кинетических уравнений Шахова [7], которые для многих задач значительно точнее круковской модели. Например, в третьем приближении аппроксимирующее уравнение имеет вид

dfjdt = ч (/+ —/), (7)

где /+ = /о^1 +(1 + Рг) (*2^7^----------1-)] ’ Рг — числ0 Прандтля,

qt—-поток тепла, с( — тепловая скорость.

В этом случае скорости молекул после столкновения распределены в соответствии с функцией /+.

2. Приложение метода к решению осесимметричных задач.

При решении задач для плоских течений траектории частиц есте-

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

d-f_ д/ | i. df | j. df / о\

dt ~ dt Kx дх + КУ dy ’ W

и координаты частицы в процессе свободного пролета за время Д*, где

а 1 ,/ Ах Ду

Д^ = ш1п'

определяются по формулам

16*1 1 16»

-)•

Хп+1 — хп Д£, Уп+1 =уп + гуЫ.

(9)

При решении осесимметричных задач траектории частиц строятся в пространстве и, вообще говоря, в формуле (8) необхо-

* д/

димо учитывать член •

Для того, чтобы использовать симметрию течения по азимутальному углу, в работе [3] был предложен следующий прием. В поле течения вводится цилиндрическая система координат, в которой ось х совпадает с осью симметрии. Поле течения разбивается на ячейки по х, г и <р. И поскольку во всех ячейках по течение одинаковое, траектории частиц строятся только в одной ячейке по ср. Таким образом, из цилиндра, который представляет собой поле течения, вырезается долька с углом раствора Д<р. Траектории частиц берут начало и оканчиваются на границах течения — торцах и верхней части дольки, а на ее гранях задается закон зеркального отражения частиц. Такой прием хоть и не избавляет

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

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

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

*

<*/_ д/ , е д/ , е д/ , $ д/ М, д/ /1ЛЧ

■й-“ ^ + ~д^ + ~г~ж;—г-щ •

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

г, но в отличие от прямолинейных траекторий (9) в нашем случае

траектории будут искривляться под действием сил

Е2 Е £

--т*. (П>

И хотя координаты смещения частиц за время М будут определяться по формулам, аналогичным (9) —

хп+\ = хп + ЪХМ\ гп+1 = гп 4- %гпМ,

скорости частиц за?это же время изменятся под действием сил (11)

(12)

, £ (£г Ол Д .

---- ?(рп --------------- --------- •

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

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

В случае плоских течений погрешность вычислений на п-м шаге дается формулой (5) .

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

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

Во втором случае для оценки погрешности условия (5) уже недостаточно, поскольку бесстолкновительное движение частиц совершается в поле сил (11) и сопровождается постоянным изменением скоростей (12). Таким образом, необходимо дополнительное условие на точность вычисления скоростей. Погрешность вычисления скорости за п шагов для уравнения

8Фп<8Фо + лд^О(д^ д*. ДУ)-В случае осесимметричных течений

8фп<;8фо + пМО(А(, Ах, Дг, Д?)

(13)

О^Тж|~яДЮ(Д*, Ах’ Аг’ А^'

определяется по формуле [4]:

0(Д*)2 + (1+ДШ)^Го,

(14)

где

Раскладывая (14) по малым величинам М, получим

Поскольку здесь число п невелико, в практических расчетах не превосходит десяти (столько ячеек проходит частица между столкновениями), то дополнительным условием при оценке погрешности метода можно считать условие

Таким образом, Д£ при больших г определяется из условия устойчивости (4), в случае же, когда траектория частицы проходит вблизи оси симметрии, М выбирается с учетом условия (15). Интересно отметить, что в случае, когда моделирование течение рассматривается в ячейке с углом раствора Д<р и зеркально отражающими гранями, выбор М происходит автоматически.

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

На боковой и задней поверхности области в качестве граничного ставится условие отсутствия градиента функции распределения.

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

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

где ^оо—средняя длина свободного пробега молекул между столкновениями на бесконечности, р0—коэффициент вязкости при тем-

Щг<сі-

(15)

т

Кп = Хоо/£, Ке„ =

ГО

пературе торможения, /. — характерный размер тела и безразмерная макроскопическая скорость 5, аналог числа М

2 кТ

Для получения результата с погрешностью в несколько процентов, как видно из формулы (13), количество траекторий должно быть порядка десятка тысяч, что и подтверждается расчетами работы [3], где было показано, что для установления течения достаточно 30 тыс. траекторий. Для вычисления с той же точностью интегральных характеристик, таких как сопротивление, количество частиц может быть на порядок меньше.

К'1=1,м~^Т,ти/=тв

5=5 ,Нп=0^,8= 18\ (1~Т

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

В качестве характерных примеров картины полей течений на фиг. 1 и 2 представлены изолинии плотности в окрестности сферы и конуса. На фиг. 3 показано изменение поля плотности на поверхности конуса в зависимости от числа Кп, а на фиг. 4 и 5 представлен характер изменения сил сопротивления для сферы и конуса по мере изменения разреженности среды. Результаты расчетов по данной методике сравниваются с результатами расчетов других авторов и экспериментов, и это сравнение показывает удовлетворительное согласие с имеющимися в литературе результатами. Звездочками на графиках обозначены результаты данной работы. На фиг. 4 значения коэффициента сопротивления сферы с температурой поверхности, примерно равной температуре на бесконечности при 5=10, сравниваются с экспериментальными

данными (штрихованная область), взятыми из работы [9]. Эксперименты проводились при гиперзвуковых скоростях. На фиг. 5 результаты расчетов сх конуса при нулевом угле атаки сравниваются с результатами расчетов по приближенной теории [10] (пунктир). Сплошной прямой на фиг. 5 показано значение невязкого предела.

2,4

2,0

¥

гг

0,8

2,0

1,0

’ 10'1 10° 10 10г 103 Яеа

Фиг. 4

. ?;.V* в — > -10е г— — *

/ / k S’ s

ю~* 1СГ3 ю* ж1 ІР W' /171 =

Фиг. 5

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

ЛИТЕРАТУРА

1. Ларина И. Н. Обтекание сферы разреженным газом. ПММ, т. 33, вып. 5, 1965.

2. V о g е n i t е F. W., Bird G. А., BroadwellJ. Е., Rungal-d i е г Н. Theoretical and experimental study of rarefied supersonic flows abent several simple shapes. .AIAA J.“, N 12, 1968.

3. Хлопков Ю. И. Статистический метод решения приближенного кинетического уравнения. „Ученые записки ЦАГИ*, т. 4,

№ 4, 1973.

4. Дьяченко В. Ф. Основные понятия вычислительной математики. М., .Наука”, 1972.

5. Хэвиленд Дж. К. Решение двух задач о молекулярном течении методом Монте-Карло. В сб. .Вычислительные методы в динамике разреженных газов", М., „Мир”, 1969.

6. Н о 1 w а у L. Н. New statistical models in kinetic theory: methods of construction. ,Phys. Fluids”, vol. 3, N 3, 1966.

7. Шахов Б. М. Метод аппроксимации кинетического уравнения Больцмана. В сб. .Численные методы в теории разреженных газов”, ВЦ АН СССР, 1969.

8. Шахов Б. М. Метод исследования движений разреженного газа. М., .Наука”, 1974.

9. KussoyM. I., Horstman С. С. Cone drag in rarefied hypersonic flow, AIAA J., vol. 8, N 2, 1970.

10. Ortloff C. R. Hypersonic low density transitional regime flow over conical vehicles. „The Journal of the Astronaut Scienses”, vol. XVI,

N 2, 1969.

Рукопись поступала 21 / VI 1977 г.