CC BY
134
Системный анализ
41
УДК 533.6
Ю. В. ЯЦКЕВИЧ, Ю. В. КОЖЕДУБ, Белорусский национальный технический университет, Минский государственный высший авиационный колледж
КОМПЬЮТЕРНАЯ МОДЕЛЬ ОБТЕКАНИЯ ЛЕТАТЕЛЬНОГО АППАРАТА ВОЗДУШНЫМИ ПОТОКАМИ И РАСЧЕТ ЕГО АЭРОДИНАМИЧЕСКИХ ХАРАКТЕРИСТИК
В работе описывается оригинальная полуэмпирическая модель турбулентного течения воздуха в пограничном слое на поверхности летательного аппарата. На основе этой модели разработана компьютерная система, позволяющая моделировать аэродинамические процессы вокруг летательного аппарата и по результатам моделирования рассчитывать его аэродинамические характеристики. Эти характеристики используют авиаконструкторы для оценки качества и свойств летательного аппарата.
The original semiempirical model of turbulent flow of air in an interface of an aircraft surface is described in this work. Based on this model the computer system allowing to simulate aerodynamic processes around an aircraft and to calculate it’s aerodynamic characteristics was developed. These characteristics are used by aircraft designers for an assessment of quality and properties of an aircraft.
Современный уровень проектирования летательных аппаратов подразумевает использование тех или иных методов компьютерного моделирования аэродинамической среды, обтекающей За?-модель планера при его движении.
Поиск приемлемых для практики форм математического описания турбулентных течений или моделей турбулентности идет уже на протяжении более 100 лет (начиная с классических работ О. Рейнольдса). Это объясняется исключительной сложностью турбулентности как физического явления.
Турбулентные течения, согласно современным представлениям, подчиняются классическим уравнениям Навье-Стокса, и в этом смысле проблема может считаться давно решенной. Однако, несмотря на фантастический прогресс вычислительной техники, наблюдаемый в последние десятилетия, ее возможности все еще недостаточны для решения этих уравнений при высо-
ких разрешениях расчетной сетки (малых шагах дискретизации пространства), представляющих практический интерес. Такая ситуация по самым оптимистичным прогнозам останется вплоть до второй половины, а то и до конца XXI века [1]. В связи с этим, как и ранее, ключевым вопросом в рассматриваемой области является поиск приемлемого компромисса между физической адекватностью модели и приемлемым для практического применения уровнем ее сложности.
Данная работа посвящена оригинальному программному обеспечению («АэроПро-1»), построенному на основе разработанной авторами относительно простой полуэмпирической «подсеточной» модели турбулентного течения на основе уравнения Рейнольдса. То есть, расчет распределения скорости u(x, y, z, t) по трехмерному пространству моделируемой воздушной среды производился на основе уравнения Рейнольдса [2]:
дих дих дих дих
+ их ±. + иу ±. + uz ±.
дл х дх y ду z Dz
ди y ди y ди y ди y
+ их ■ + и,. - + и7
д( х дх y ду z Dz
ди_ ди_ ди7 Duz
+ ит + иУ + и7 ^ ■
д1 дх У ду z Dz
1 dP и
-----+ _•
р дх р
(я 2 я2
д и^ д и^ д
дх2
+
дy
+
и
дz 2
ди'х и'у
дy
1 дР и
-----+—-
pdy р
1 дР и
----+ —
р дz р
( У2 ~,2 ~2
д и,, д и,, д и
дх2
+
V
( У2
ду2
+
дz 2
ди' y и 'х
ди'х и 'z
дz
ди' y и 'z
дх
дz
д 2 и
дх
~,2 ~0. ^
z + д uz + д uz
ду
2
дz
2
ди' y и 'z
ди'х и'z
дy
дх
(1)
1, 2015
СИСТЕМНЫЙ АНАЛИЗ И ПРИКЛАДНАЯ ИНФОРМАТИКА
42
Системный анализ
где t - время; x, у, z - декартовы координаты (ось Y направлена вертикально); их, иу, uz -проекции скорости воздуха на координатные оси X, Y, Z соответственно; P - давление воздуха; р - плотность воздуха; h - динамическая вязкость воздуха; u’x, и’у, u’z - пульсации соответствующих проекций скорости воздуха внутри сеточного элемента.
Расчет распределения плотности р(х, y, z, t) на основании трехмерного поля скоростей u(x, y, z, t) производился по уравнению неразрывности:
5р 5р 5р 5р
—+ их — + и..— + uz — dt dx y dy dz
= р
(2)
Пересчет распределения плотности р(х, у, z, t) в поле давлений Р(х, у, z, t) производился по уравнению Менделеева-Клапейрона:
p=-р RT, (3)
M
где Т - температура воздуха; М = 0,029 кг/моль -молярная масса воздуха; R = 8,3144 Дж/ (моль-К) - универсальная газовая постоянная.
Данная система уравнений (1-3) является не замкнутой. Для их замыкания нужно определить 6 компонент турбулентных (Рейноль-дсовых) напряжений и \ и'}- .
Для определения Рейнольдсовых напряжений используется понятие пограничного слоя. Пограничный слой - тонкий по сравнению с характерным линейным размером тела слой жидкости или газа, прилегающий к твердой поверхности, в котором градиенты газодинамических переменных в нормальном к стенке направлении столь велики, что инерционные силы и силы трения имеют здесь один и тот же порядок. Течение в этом слое при больших числах Рейнольдса становится турбулентным. Понятие пограничного слоя [3] для анализа движения жидкости при больших числах Рейнольдса было предложено Л. Прандтлем (1904).
Разработанная модель основывается на эмпирических наблюдениях за поведением турбулентного слоя, смоделированного методом прямого численного моделирования (Direct Number Simulation - DNS) для малого участка поверхности (1x5 см). При этом использовал-
ся шаг сетки от 0,05 до 0,2 мм. DNS-модель построена на основе уравнений Навье-Стокса, уравнения неразрывности и уравнения Менделеева-Клапейрона [4]. Эта модель полностью свободна от эмпиризма (в данной статье не рассматривается). Она позволила выявить основные закономерности возникновения и развития турбулентного слоя в зависимости от шероховатости поверхности, скорости и угла набегающего ламинарного потока по отношению к обтекаемой поверхности (локального угла атаки).
Упрощенно предполагается, что пограничный слой до момента отрыва состоит из движущихся цилиндрических вихревых областей (вихрей), диаметр 5W которых равен толщине слоя 5. Линейная скорость воздуха прилегающего к твердой поверхности равна 0. Воздух на внешнем крае пограничного слоя (внешняя часть вихрей) имеет скорость набегающего ламинарного потока и. Таким образом, центр вращения вихря (и сам вихрь) перемещается со скоростью = u/2 в направлении воздушного потока.
Исходя из этого представления, Рейноль-дсовы напряжения зависят от толщины пограничного слоя (диаметра вихрей) и скорости ламинарного потока (усредненной скорости) приблизительно так:
1 8 3 Ay
1 Ay
3 8х
>2
У
л2
при 8 < Ау,
при 8 > Ау,
(4)
где u’x - пульсация скорости воздуха вдоль обтекаемой поверхности; и’у - пульсация скорости воздуха вдоль нормали к обтекаемой поверхности; и 'x и'у - Рейнольдсово напряжение; 5 - толщина пограничного слоя; ux - проекция скорости ламинарного потока на обтекаемую поверхность; Ау - высота сеточного элемента.
Как показало DNS-моделирование, турбулентный вихрь появляется на поверхности при наличии неровности и еще двух условиях:
1) набегающий поток направлен от поверхности, то есть локальный угол атаки внутри сеточного элемента не отрицательный (ai > 0);
2) скорость набегающего воздушного потока достаточно велика (и > и0).
СИСТЕМНЫЙ АНАЛИЗ И ПРИКЛАДНАЯ ИНФОРМАТИКА
1, 2015
Системный анализ
43
Рис. 1. Панель задания входных параметров при создании новой модели
Значение и0 зависит от шероховатости поверхности. В данной работе принимается и0 = 7 м/с, что соответствует наличию на поверхности неровностей высотой 0,4 мм.
По мере перемещения вихря вдоль набегающего потока наблюдается рост его геометрических размеров в результате вовлечения во вращательное движение окружающих масс воздуха. Этот процесс характеризуется скоростью роста толщины турбулентного пограничного слоя (Kg = A5/At). По результатам численных экспериментов с различными скоростями набегающего потока была установлена средняя величина скорости роста толщины турбулентного потока - Kg « 2 м/с.
Кроме роста размеров вихрей на геометрическую конфигурацию турбулентного слоя оказывает влияние усредненная скорость течения воздуха. В целом изменение толщины турбулентного слоя в разработанной модели описывается уравнением:
dg dg dg dg
— + ux— + uy— + uz— = Kg . (5)
dt dx y dy dz
Таким образом, уравнения (1-5) представляют собой замкнутую систему, описывающую динамику воздушных потоков в пространстве вокруг летательного аппарата. В нее также необходимо включить условия образо-
вания турбулентного слоя (aj > 0 и u > Uo), то есть в компьютерной модели анализируется угол наклона набегающего потока к обтекаемой поверхности в каждом граничном конечно-разностном элементе.
Эта математическая модель была положена в основу компьютерной модели, в которой моделируемое пространство разбивается на прямоугольные элементы (прямоугольную сетку). Элементы, пересекающие поверхность твердого объекта, имеют трапециевидную форму (конечно-объемное представление на границах твердого объекта).
Кроме описания геометрии объекта в формате «*.stl>> входными данными для моделирования являются: скорость набегающего потока, угол атаки (угол между направлением воздушного потока и продольной осью самолета в вертикальной плоскости), угол скольжения (угол между направлением воздушного потока и продольной осью самолета в горизонтальной плоскости), размеры воздушного пространства с разных сторон относительно летательного аппарата, а также - характеристики атмосферы (рис. 1).
На рис. 2, 3 представлены результаты моделирования в виде полей скоростей и давлений воздуха по сечениям моделируемого пространства, а также в проекциях на поверхность самолета.
1, 2015
СИСТЕМНЫЙ АНАЛИЗ И ПРИКЛАДНАЯ ИНФОРМАТИКА
44
Системный анализ
Рис. 2. Интерфейс программы при отображении распределения скоростей воздуха по сечениям моделируемого пространства и по поверхности самолета при угле атаки набегающего потока a = 0°
(£*-5.СТ5м_ iZ-J.«®83 M jttami™:-32.S1S3 jBm-»: 0.11D7«c
Рис. 3. Интерфейс программа: при отображении распределения давлений воздуха по сечениям моделируемого пространства и по поверхности самолета при угле атаки набегающего потока a = 0°
На основе полученных трехмерных полей давлений рассчитываются основные аэродинамические характеристики самолета [5] (лобовое сопротивление и подъемная сила, моменты тангажа, рыскания и крена).
Аэродинамические характеристики вычислялись следующим образом. Вначале вычислялись силы давления, действующие с разных сторон на самолет, путем суммирования сил давления каждого конечно-объемного
СИСТЕМНЫЙ АНАЛИЗ И ПРИКЛАДНАЯ ИНФОРМАТИКА
1, 2015
Системный анализ
45
элемента, пересекающего поверхность объекта:
Fn= Е (P(x,y,z)-P0)AyAz , (6)
(x,y,z) - воздух (x+Ax,y,z) - тело
F3= Е {.P{x,y,z)-P0)AyAz, (7)
(x,y,z) - воздух (х-Аx,y,z) - тело
Fe= Е (P(x,y,z)-P0)AxAz, (8)
(x,y,z) - воздух (x9y—Ay,z) — тело
FH= Е (P(x,y,z)-P0)AxAz , (9)
(x,y,z) - воздух (x,y+Ay,z) - тело
где Fn - сила давления, действующая на объект спереди, F3 - сзади, FB - сверху, FH - снизу, P(x, y, z) - давление в конечно-разностном элементе с координатами (x, y, z); P0 = 105005 Па -нормальное давление воздуха; Ax, Ay, Az - шаг конечно-разностной сетки.
Затем, на основании полученных сил, действующих с разных сторон на объект, вычислялись лобовое сопротивление:
X = Fn~F3, (10)
подъемная сила:
7 = FH-FB, (11)
коэффициент лобового сопротивления:
С - 2Х ^ ^нР0и»2’
коэффициент подъемной силы: 27
Су =-------2,
^нР0исю
(12)
(13)
где Ан - площадь несущей поверхности (вертикальной проекции самолета), Р0 - нормальная плотность воздуха (на бесконечном удалении от самолета); - скорость воздуха на бесконечном удалении от самолета.
Для вычисления вращающих моментов нужно более подробное описание системы координат: ось Х направлена вдоль оси самолета от носовой части к хвостовому оперению, ось Y направлена вертикально (перпендикулярно плоскости крыльев), ось Z направлена вдоль крыльев. Тогда крен соответствует вращению самолета вокруг оси X, рыскание - вращение самолета вокруг оси Y, тангаж - вращение самолета вокруг оси Z.
В разработанной программе вращающие моменты вычисляются относительно осей вращения, проходящих через центр тяжести самолета, координаты которого (хс, yc, zc) задаются в программе перед началом моделирования. В процессе расчета просматриваются все пары соседних конечноразностных элементов, в которых один элемент принадлежит твердому телу (летательному аппарату), а второй - воздушному пространству. Произведение давления в пограничном воздушном элементе на его расстояние до центра тяжести и площадь грани дает вращающий момент. А сумма всех элементарных вращающих моментов дает полный вращающий момент летательного аппарата:
Мх = X Р(Х У, z)(z - zc ) AxAz -
(x,y,z) - воздух (x,y+Ay,z) - тело
- X P(x у, z)(z - zc ) AxAz +
(x,y,z) - воздух (х,у-Ау,г)-тело
+ E P(x> 7, z)(y - ус )Д*Ду -
(x,y,z) - воздух (x,y,z+Az) - тело
- E P(X У, z)(y - yc )AxAy,
(x,y,z) - воздух (x,y,z+Az) - тело
My= X P(x,y,z){z-zc)AyAz-
(x,y,z) - воздух (x+Ax,y,z) - тело
- E p(x> у, z)(z - zc )ау^+
(x,y,z) - воздух (x+Ax,y,z) - тело
+ X P(F y, z)(x - xc ) AxAу -
(x,y,z) - воздух (x,y,z+Az) - тело
- E p(x’ У> z)(x ~ xc )AxAy,
(x,y,z) - воздух (x,y,z+Az) - тело
mz = E p(x> y, z)(y - yc)Aу*2 -
(x,y,z) - воздух (x+Ax,y,z) - тело
- E p(x y, z)(y - yc )ау^+
(x,y,z) - воздух (x+Ax,y,z) - тело
+ X p(x> y>z)(z~zc )AxAz -
(x,y,z) - воздух (x,y+Ay,z) - тело
(16)
- X P(x,y,z)(z-zc)AxAz,
(xfy,z) - воздух (x,y+Ay9z) — тело
где Mx, My Mz - моменты крена, рыскания и тангажа соответственно; P(x, y, z) - давление воздуха в конечно-разностном элементе c коорди-
1,2015
СИСТЕМНЫЙАНАЛИЗ И ПРИКЛАДНАЯ ИНФОРМАТИКА
46
Системный анализ
^еродчнамчмеские характеристмси
Угол атаки:|0
Спереди
, скольжения:|0 Сзади Сверху
Скорость потока (м/с): 33
Снизу Справа
Слева
(п) ( з ) ( Б ) (н) (пр) (л)
Давление (Р): 117.828 -37.7G5 -166.008 -59.2153 -48.4413 -50.0187
Площадь (S): 3.25255 3.25255 19.6746 19.8746 8.80234 8.80234
Сила (F=PKS): 383.246 -122.832 -3266.17 -1165.04 -479.683 -495.302
Вдоль
потока (Хтр): 62.0781
Вертикально потоку (Утр): -0.445944
Поперек
потока (Zip): 0.00571129
[Лобовое сопротивление (X=Fn+F3+XTp): 508.078 Подъемная сила (Y=Fb+Fh+Ytp): 2101.13 Боковая сила (Z=Fnp+Fn+Zrp): -15.6142
|&®ФФ. лобового [сопротивления (Сх):0.03856 КЙ®ФФ. подъемной [ силы (Су): 0.1801 Коз ФФ. боковой силы (Сг): -0.00118
[Аэродинамическое качество (Y/X): 4.152 |
Коорд.центра масс - X: |-5.5 Y: |-2.В5 Z: |4.837 Размах крыла(м):|4.834
Моменты
давления:
трения:
Сумм, момент: Коэффициент момента:
Крен Рыскание
(Мрх): 2G.2444 (Мру): -56.3368
[Мук]: -0.01214: (Mvy): 0.596255
[Мх): 26.2323 [Mil): -55.7505
[тх): 0.000413 (ту): -0.0003Л
Т ангаж
(Mpz)
(Mvz)
(мГ
-507.118
-5.09221
-512.21
Форм.Файл |
(mz): -0.00807'
OK
Рис. 4. Панель результатов расчета аэродинамических характеристик программой «АэроПро-1»
натами (x, y, z); Ax, Ay, Az - шаг конечно-разностной сетки.
На основании полученных моментов вычислялись соответствующие коэффициенты:
Мх
тх = 2, ,
*^нРоиоо I
Му
mY=-----
-ShPoO
м7
mz =----
*^нРоиоо I
(17)
(18) (19)
где mx my mz - аэродинамические коэффициенты моментов крена, рыскания и тангажа соответственно; «Х - площадь несущей поверхности (вертикальной проекции самолета), р0 -нормальная плотность воздуха (на бесконечном удалении от самолета); - скорость воз-
духа на бесконечном удалении от самолета; l -размах крыльев.
Пример результатов расчетов аэродинамических характеристик самолета программой «АэроПро-1» представлен на рис. 4.
Точность расчетов программы оценивалась путем моделирования и сравнения основных аэродинамических характеристик крыльев типовых профилей (лобовое сопротивление и подъемная сила) с известными данными [6] по результатам продувки этих крыльев в аэродинамических трубах. Это сравнение показало, что разработанная компьютерная модель количественно и качественно хорошо имитирует реальное поведение крыльев на углах атаки до 120. Погрешность определения лобового сопротивления не превышает 13%, а погрешность определения подъемной силы не превышает 5%. Однако, на больших углах атаки (>120), соответствующих срыву воздушного потока с обтекаемой поверхности, наблюдается значительное несоответствие с натурными экспериментами.
Литература
1. Гарбарук А. В. Моделирование турбулентности в расчетах сложных течений: учебное пособие / А. В. Гарбарук, М. X. Стрелец, М. Л. Шур - СПб: Изд-во Политехи. ун-та, 2012. - 88 с.
2. Андерсон Д., Таннехилл Дж., Плетчер Р Вычислительная гидромеханика и теплообмен/ В 2х т., т. 2: [Пер. с англ.]. - М.: Мир., 1990.
3. Шлихтинг Г. Теория пограничного слоя, пер. с нем., М.: Наука, 1974. - 711 с.
4. Белов И. А., Исаев С. А. Моделирование турбулентных течений: Учебное пособие / Балт. гос. техн. ун-т. СПб., 2001. 108 с.
5. Микеладзе В. Г., Титов В. М. Основные геометрические и аэродинамические характеристики самолетов и ракет: Справочник. - М.: Машиностроение, 1982. - 149 с.
СИСТЕМНЫЙ АНАЛИЗ И ПРИКЛАДНАЯ ИНФОРМАТИКА
1, 2015
