_______УЧЕНЫЕ ЗАПИСКИ Ц АГ И
Т о м IX 197 8
№ 5
УДК 533.6.011.8
ПРОСТРАНСТВЕННОЕ ОБТЕКАНИЕ ПЛАСТИНЫ ГИПЕРЗВУКОВЫМ ПОТОКОМ РАЗРЕЖЕННОГО ГАЗА
А. И. Ерофеев
Приводятся результаты расчета методом Монте-Карло аэродинамических характеристик и поля течения около прямоугольной и треугольной пластин при значениях чисел = 10 и 20 и 1*е0<;30, углах атаки 0 и 15°. Расчеты проводились для одноатомного газа. Сечение столкновения молекул принималось постоянным (молекулы— твердые шары). Результаты сопоставляются с расчетными данными для случая пластины бесконечного размаха.
1. Метод расчета. Решение уравнения Больцмана проводилось одним из вариантов метода Монте-Карло -— методом прямого моделирования. Поскольку описание этого метода подробно дано в ряде работ, например [1—3], здесь ограничимся кратким изложением метода и особенностей программы расчета. В методе прямого моделирования прослеживается движение ансамбля молекул, заменяющих реальный газ. Около обтекаемого тела выделяется некоторая область, которая разбивается на ячейки (фиг. 1) с линейным размером Л, меньшим местной длины свободного пробега молекул. В каждую ячейку в начальный момент времени помещается определенное число молекул, скорости которых соответствуют начальному виду функции распределения (произвольной).
Задача решается методом установления, процессы движения молекул и столкновения между ними рассматриваются последовательно в течение шага по времени М. Столкновения между молекулами описываются статистически, причем сталкиваются только молекулы, принадлежащие одной геометрической ячейке. За шаг М с границ области проводится „вбрасывание" определенного числа молекул в соответствии с граничной функцией распределения. При своем движении в расчетной области молекулы могут столкнуться с телом—-в этом случае записывается информация о принесенных на поверхность тела импульсе и энергии. При суммировании этой информации по большому числу столкновений
молекул с поверхностью (обычно число столкновений составляет несколько тысяч) определяются суммарные и локальные аэродинамические характеристики.
В различных зонах течения длина свободного пробега молекул может существенно различаться. Поскольку в методе прямого моделирования размер ячеек А должен быть меньше местной длины свободного пробега, то при равномерном разбиении расчетной
области размер ячеек должен устанавливаться по наименьшей длине свободного пробега Тогда в областях течения, в которых
разбиение оказывается излишне мелким. В связи с этим в программе была предусмотрена возможность неравномерного разбиения области течения (см. фиг. 1): в зоне 1 около пластины выбиралась сетка с шагом ки. причем хф Л, у ф А, г, в остальной области — зона 2—разбиение проводилось с шагом А, большим и одинаковым во всех направлениях. Такое простое разбиение области было выбрано для того, чтобы избежать дополнительных затрат расчетного времени при определении принадлежности молекул данной ячейке, которые будут иметь место при более сложном разбиении.
Введение неравномерной сетки позволяет повысить точность расчетов в областях течения с большими градиентами. При обтекании пластины под малыми углами атаки такими областями являются зона вблизи поверхности (особенно в случае холодной стенки) и зона формирования ударной волны, хотя, конечно, в разреженном газе фронт скачка уплотнения сильно размыт. Поэтому было введено более мелкое разбиение для области вблизи поверхности пластины.
Однако уменьшение объемов ячеек приводит к уменьшению числа молекул, находящихся в них на каждом шаге по времени, что может привести к снижению, точности расчетов. Поэтому было предусмотрено введение различных статистических весов для молекул в указанных зонах течения (метод „расщепления" и „рулетки").
Полагается, что при пролете молекулы из зоны 2 (размер ячейки к) в зону 1 (размер ячейки Л,) происходит расщепление траектории на две, и каждая молекула приобретает статистический вес /э1= 1 /2; при пролете в обратном направлении каждая вторая траектория обрывается, а оставшаяся молекула приобретает статистический вес Р= 1. Возможность неравномерного разбиения области с введением статистических весов была проверена в серии контрольных расчетов на равномерной и неравномерной сетке для одинаковых условий обтекания (см. п. 2). Отметим, что применение такой методики позволило продвинуться в расчетах в сторону больших значений числа Ие0 (меньших значений числа Кнудсена).
На фиг. 1 приведена картина расчетной области и геометрической конфигурации пластины. Граничные условия ставятся следующим образом: на границах вверх по потоку — на плоскости х = 0 и при афО на плоскости у = 0 — функция распределения предполагается невозмущенной, на плоскости симметрии 2 = 0 и при угле атаки а = 0 на плоскости у = 0, не занятой пластиной, ставится условие зеркального отражения молекул. На остальных границах области принимается условие отсутствия градиента функции распределения. Функция распределения молекул, отраженных от поверхности пластины, предполагается максвелловской с температурой, равной температуре поверхности Тт. В общем случае пластина имеет трапециевидную форму в плане,, причем точки 3 и 4 могут задаваться произвольно на линиях х = х{ и х = х2.
2. Результаты расчетов. Расчет обтекания прямоугольных и треугольных плоских пластин проводился для одноатомного газа (молекулы — твердые шары) для чисел Моо=10 и 20 и двух значений температурного фактора іш = Тт1Т0 — 1 и 0,1 (7'0 — температура торможения). Максимальное число Ие0 не превышало 30: Ие0 — Роо ^СО ^ гг
= —,т Л , р.», 1Л» — плотность и скорость газа в невозмущенном
Р- V ‘ о)
потоке, Ь — длина пластины вдоль по потоку, р(Т0) — коэффициент вязкости, вычисленный по температуре торможения.
Важным параметром в методе прямого моделирования является количество молекул, находящихся в геометрической ячейке. Контрольные расчеты при а = 0 показали, что интегральные характеристики слабо меняются при изменении числа молекул в невозмущенной области течения М* от 2 до 10: для Не0 = 3,76, *„,= 1, 1г/Ь = 2 (1г — размах пластины вдоль оси г), вариации коэффициентов сопротивления сх и нормальной силы сп не превышали 5%, а при изменении УУ» от 5 до 10 — менее 2%. Однако с изменением числа молекул менялась частота столкновений. Так, при Л/оо=Ю расчетная частота столкновений в ячейках в невозмущенной области течения была близка к теоретическому значению, при А^оо = 5 различие в расчетной и теоретической значениях частот столкновений составляло примерно 10%, а при = 2 — примерно 20%. Поэтому при проведении расчетов количество молекул выбиралось так, чтобы в ячейках в возмущенной области около пластины Их число было примерно равно 5.
Контрольные расчеты для выяснения возможности введения неравномерной сетки и различных статистических весов проводились при обтекании треугольной пластины с углом при вершине 26 = 34° для следующих значений параметров: 1) М*, = 10, Нео=20, £го=1; 2) Ке0=Ю, Моо^Ю, ^ = 0,1. В обоих случаях различия в сх и сп для равномерной и неравномерной сеток не превышало
4%. Зона 1 при этом простиралась на расстояние от пластины, равное двум-трем длинам свободного пробега молекул. Введение неравномерной сетки дало возможность продвинуться в расчетах в сторону больших чисел Ие0. Это можно проиллюстрировать на примере обтекания треугольной пластины 2в = 53° при а=15°, Моо = 20, <то=1, Не0 = 20. При расчете на равномерной сетке с шагом Л величина сх была равна 0,41. Если теперь „затрубить* сетку до размера Л* = 1,67Л, то величина сх изменяется на 15%. Введение же неравномерной сетки с х = 0,75 Л*, Ь1х = А, у = 0,5 А* и статистических весов уменьшает ошибку в сх до 4%. Поэтому расчеты при максимальных числах Яе0 проводились на неравномерной сетке с применением методики „расщепления11 и „рулетки", в этом случае расчет на равномерной сетке с мелким шагом невозможен из-за ограниченности оперативной памяти БЭСМ-6, а расчет на равномерной сетке с крупным шагом приводит к большим ошибкам при вычислении аэродинамических коэффициентов.
Статистическая погрешность вычислений аэродинамических коэффициентов оценивалась по флуктуациям величин в установившемся режиме обтекания. Осредненные по шагам значения сх и су выводились на печать через каждые 100 шагов (полное число шагов в квазистационарном режиме было порядка 1000). При а = 0 максимальные флуктуации коэффициента сопротивления не превышали +2%, при а =15 и 30° флуктуации сх и су были менее +1,5% при
Таблица 1
При М00= 11, <„,= ], Яе0 = 3,76
ЬП 1 2 4 со
Сх 0,063 0,063 0,065 0,068
Сп 0,05 0,054 0,058 0.С67
Таблица 2
При Мю= 10, tw = I, Иео— Ю
ць 1 2 оо
сх 0,074 0,082 0,0865
сп 0,073 0,086 0,098
всех расчетных числах Ие,,. При этом, конечно, расчеты при больших числах Ие^ требовали больших затрат времени из-за увеличения области течения, т. е. из-за увеличения числа молекул, моделирующих реальный газ. Что касается статистической погрешности при вычислении поля плотности (и других параметров газа),^го она может быть оценена по порядку величины как Ьп — 0(1/|Л/У), где ./V—полное число частиц в ячейке, просуммированное по всем шагам по времени к в квазистационарном режиме течения (ДГ=. пЫоок). Учитывая приведенные выше значения величин к, полное число частиц в ячейках было порядка нескольких тысяч, т. е. погрешность в вычислении плотности можно оценить в несколько процентов.
В табл. 1 и 2 приведены результаты расчета аэродинамических характеристик прямоугольной пластины, показывающие влияние относительного удлинения.
Как следует из приведенных данных, при небольших числах Ке0 размах пластины становится существенным параметром, особенно сильно влияющим на величину давления.
На фиг. 2 и 3 приведены результаты расчета коэффициентов сопротивления с* и нормальной силы с*п, отнесенных к своим свободномолекулярным значениям, при ^ = 0, Мсо=Ю, Ьт=\ (фиг. 2) и ^ = 0,1 (фиг. 3) для треугольной пластины (2 0 = 34°) — сплошные кривые, и для пластины бесконечного размаха — пунктирные кривые. Уменьшение аэродинамических коэффициентов связано с уменьшением возмущения поля течения в пространственной задаче по сравнению со случаем, плоского течения. Возможно, что длина
пластины, по которой рассчитывается число Рейнольдса, не яв-
Сх
1,50
1,25
1,00
3,0
2.5 2,0
1.5
/
/ /
/ / /
/ / /
/ / /
S" / / J
.
\ \ \
/ /
/ / /
3 10 Ре0
Фиг. 2
J 10 Фиг. 3
А?0
ляется характерным размером в этом случае. Но из фиг. 2 и 3 следует, во-первых, что такой характерный размер зависит от температурного фактора, а, во-вторых, вид зависимости c,,.(Re0) показывает, что при tw = 1 согласования результатов расчетов для плоской и пространственной задач не может быть достигнуто над* лежащим выбором характерного размера.
На фиг. 4 приведены поля плотности в некоторых сечениях х — const вблизи плоскости симметрии z — 0 при обтекании прямоугольной пластины: Мю = 10, Re0 = 10, IJL=\ (кривая 1), IJL = 2 (кривая 2), пластина бесконечного размаха (кривая 3). Видно, что с уменьшением отношения IJL возмущение поля плотности уменьшается. На фиг. 5 приведены поля плотности при обтекании прямоугольной пластины (Мею =10, Re0=10, IJL = 2, = 1) в некото-
рых сечениях x/L, z'L, где z — расстояние от плоскости симметрии
б — Ученые записки № 5
81
—%— * -----\
Фиг. 4
°1 1,1 1 1,1 0,3 1 1,1
—•— г/С~0,083
Фиг. 5
2 = 0, иллюстрирующие эффект пространственности течения — при приближении к боковой кромке (увеличение г) поле ПЛОТНОСТИ возмущено меньше, чем на оси симметрии.
Аэродинамические характеристики треугольной пластины с углом при вершине 26 = 53° при обтекании под углом атаки а =15° приведены на фиг. 6 (Моо = 20)— сплошные линии. Здесь же для сравнения приведены результаты расчетов для пластины бесконечного размаха (пунктирные линии). По сравнению с обтеканием под нулевым углом атаки заметно уменьшение различий в коэффициенте сопротивления для плоской и пространственной задач. На фиг. 3 и 6 звездочками отмечены результаты расчета обтекания треугольной пластины в случаях, когда навстречу потоку обращена не вершина треугольника, а его основание. Эти
данные свидетельствуют о том, что при заданной площади и длине аэродинамические характеристики слабо зависят от формы пластины. Такое свойство течений известно в теории первых межмолекуляр-ных столкновений (теорема обратимости) [4]. Здесь же оно получено при существенно меньших числах Кнудсена, когда теория первых столкновений неприменима. Слабая зависимость коэффициента сопротивления при малых углах атаки таких тел, как пластина и конус, если за характерную площадь принимать площадь омываемой поверхности, была отмечена в работе [5].
ЛИТЕРАТУРА
1. Bird G. A. Direct simulation and the Boltzmann equation. „Phys. Fluids-, vol. 13, N 11, 1970.
2. Белоцерковский О. М., Яницкий В. Е. Статистический метод частиц в ячейках для решения задач динамики разреженного газа. Ч. 1. Основы построения метода. „Ж. вычисл. матем. и матем. физ.‘, 15, № 5, 1975.
3. Ерофеев А. И., Перепухов В. А. Расчет поперечного обтекания пластины потоком разреженного газа. „Изв. АН СССР, МЖГ“, 1976, № 4.
4. Коган М. Н. Динамика разреженного газа. М., „Наука", 1967.
5. Гусев В. Н., Крюкова С. Г. О несущих свойствах тел в переходной области при гиперзвуковых скоростях потока. .Ученые записки ЦАГИ*, т. 4, № 6, 1973.
Рукопись поступила 141X11 1977 г.