Программное обеспечение «АэроПро-1» для численного моделирования аэродинамических характеристик Текст научной статьи по специальности «Физика»

УДК 621 . 74

Поступила 08.10.2014

Ю. В. ЯЦКЕВИЧ, А. Н. ЧИЧКО, С. Г. ЛИХОУЗОВ, О. А. САЧЕК, В. Ф. СОБОЛЕВ, Т. В. МАТЮШИНЕЦ, О. И. ЧИЧКО, БНТУ

ПРОГРАММНОЕ ОБЕСПЕЧЕНИЕ «АЭРОПРО-1» ДЛЯ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ АЭРОДИНАМИЧЕСКИХ ХАРАКТЕРИСТИК

Представлено краткое описание математического аппарата для компьютерной системы «АэроПро-1», разрабатываемой для расчета аэродинамических характеристик объектов, движущихся в воздушных средах.

The short description of mathematical apparatus for the computer system AeroPro-1 developedfor calculation of aerodynamic characteristics of the objects moving in air environments is presented.

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

Цель работы - разработка специального программного обеспечения (ПО) «АэроПро-1» для моделирования аэродинамических свойств среды для решения задач проектирования геометрии планеров

Программа «АэроПро-1» разрабатывается сотрудниками Белорусского национального технического университета, имеющих опыт создания ПО «ПроЛит-1», «ПроТерм-1» и «ПроНРС-1» . В основу разрабатываемого программного обеспечения «АэроПро-1» положена математическая модель кле-точно-автоматного типа [1-4] . Для точной аппроксимации геометрии крыла разработаны клеточно-автоматные правила, позволяющие учесть аэроди-

намические особенности различных скоростных потоков, вызывающих турбулентность

Современный уровень проектирования летательных аппаратов подразумевает использование методов компьютерного моделирования аэродинамической среды, обтекающей 3d- модель планера при его движении Поиск приемлемых для практики форм математического описания турбулентных течений или моделей турбулентности идет уже на протяжении более 100 лет начиная с классических работ О . Рейнольдса . Это объясняется исключительной сложностью турбулентности как физического явления

Турбулентные течения, согласно современным представлениям, подчиняются классическим уравнениям Навье-Стокса, и в этом смысле проблема может считаться давно решенной Однако, несмотря на фантастический прогресс вычислительной техники, наблюдаемый в последние десятилетия, ее возможности все еще недостаточны для решения этих уравнений при представляющих практический интерес высоких числах Рейнольдса Даже по самым оптимистичным прогнозам будут оставаться таковыми вплоть до второй половины, а то и до конца XXI в . [5] . В связи с этим, как и ранее, ключевым вопросом в рассматриваемой области является поиск приемлемого компромисса между физической адекватностью модели и приемлемым для практического применения уровнем ее сложности . Естественно, что определение «приемлемости» существенным образом зависит от возможностей вычислительной техники, с одной стороны, и от потребностей конструктора летательного ап-

вт

пг:г г кътжпъ

4 (77), 2014-

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

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

Разработанная математическая модель клеточно-автоматного типа состоит из классической части (на основе уравнений Навье-Стокса, уравнения неразрывности, уравнения Менделеева-Клапейрона) и оригинальной части, называемой моделью пограничного слоя . Расчет распределения скорости У(х, у, z, т) по трехмерному пространству моделируемой среды производится на основе уравнения Навье-Стокса:

+дУ. ,р=о.

дг ^ дх ду дz

(2)

Пересчет распределения плотности р(х, у, z,т) в поле давлений Р(х, у, z, т) производится по уравнению Менделеева-Клапейрона:

Р = _Р.ЯТ. М

(3)

дп

дг

„ дК дУх дУх - у —^ + у —^ + у —х-

х дх у ду z дг

1 дР) { д 2у

р дх р

дк

дг

дУу

-к у

дх2 дУ„

—+уУ—

дх ду

д 2Ух

Ь-2- +

ду2

дУ

-к—

z дz

1 дР +п(Т) { д 2У

Р ду Р

дУл

дz2

12т/

у д2Уу д у

(1)

У

дг

дх2 дУ

ду2

дУ

дz2

дУ

-ух—+У у—-+У,

дх ду д,'г

1 дР +п(Т) { д 2У, Р f

д2У д2У ^

дх2 ду2

а,2

где Т - температура воздуха; М — 0,029 кг/моль -молярная масса воздуха; Я — 8,3144 Дж/(моль-К) -универсальная газовая постоянная

В конечно-разностной форме уравнения (1) и (2) имеют вид

Vx (х, у, г, г + Аг) - Vx (х, у, г, г) = А?

Р(х + Ах,у,г,г) - Р(х + Ах,у,г,г) + ц(Т(х,у,г,г)) х

— + X

р( х, у, г, г )2Ах р( х, у, г, г)

Ух (х + Ах, у, г, г) -2Ух (х, у, г, г) + У; (х -Ах, у, г, г) + Ах2

+ Ух (х, у + Ау, г, г) - 2¥х (х, у, г, г) + Ух (х, у - Ау, г, г) + Ау2

+ Ух (х, у, г + Аг, г) - 2V (х, у, г, г) + V (х, у, г - Аг, г) Ах2

-У Vx(х + Ах,у,г,г)^(х-Ах,у,г,г)

(4)

2Ах

Vx (х + Ах, у, г, г) (х - -Ах, у, г, г)

2Ау

Vx (х + Ах, у, г, г) (х- -Ах, у, г, г)

где г - время; х, у, г - декартовы координаты (ось У направлена вертикально); Ух, Уу, Уг - проекции скорости воздуха на координатные осиX, У Z соответственно; Р - давление воздуха; р - плотность воздуха; п(Т) - динамическая вязкость воздуха; Т -температура воздуха

Расчет распределения плотности р(х, у, г, т) на основании трехмерного поля скоростей У(х, у, г, т) производится по уравнению неразрывности:

г 2Аг

Vy (х, у, г, г + Аг) - Vy (х, у, г, г) = А?

= - Р(х, у + Ау, г, г) - Р(х, у - Ау, г, г) _

р(х, у, г,г )2Ау - х

р( х, у, г, г)

{Vy (х + Ах, у, г, г) - 2Vy (х, у, г, г) + Vy (х - Ах, у, г, г)

ч Ах2

V (х, у + Ау, г, г) - 2^у (х, у, г, г) + V (х, у - Ау, г, г) V—---2---+

Ау2

Уу (х, у, г + Аг, г) -2Vy (х, у, г, г) + Vy (х, у, г -Аг, г) ^

Ах2

Vy (х + Ах, у, г, г) - Vy (х - Ах, у, г, г) 2Ах

Vy(х + Ах,у,г,г)^(х-Ах,у,г,г) -2Ау

Vy (х + Ах, у, г, г) ^ (х -Ах, у, г, г) г 2Аг '

- У

-У.

У2 (X, У, t + At) - У2 (X, у, t) =

М

= - Р(X, У, 7 + М t) - Р(X, у, 7 + М-, t) + 2р( X, у, 7, t)Az

~г X

р( X, у, 7, t)

Г У7 (X + Мс, у, 7, t) - 2Уг (X, у, 7, t) + У2 (X - Ас, у, 7, t) + Х1 Ах2 +

+ V, (X, у + Ау, 7, t) - 2Vг (X, у, 7, t) + Vг (X, у - Ау, 7, t) + Ау2

Уг (х, у, 7 + А t) - 2 Vг (X, у, 7, t) + Vг (X, у, 7 - А t)

- V

-V

Ах2

Vг (X + Ах, у, 7, t) -V (X -Ах, у, 7, t) 2Ах

У2(X + Ах,у,7,t)-Уг(X-Ах,у,7,t) -2Ау

- V у (х+Aх, y, ^ ) - у (х - Aх, y, ^ t) (6)

г 2А7 '

р(X, у, 7, t + At) -р(X, у, 7, t) = Ат

Ух (X - А>, у, 7,0-Кх (х + А>, у, 7, t)

(7)

2Ах

Уу (х, у - Ау, г, г) - Уу (х, у + Ау, г, I)

2Ау

Уг (х, у, 7 - А?, г) - У2 (х, у, 7 + А г) 2А2

р( X, у, 7, t),

где Vх(х, у, 7,т), Уу(х, у, 7,т), Vz(х, у, 7,т) - проекции скорости воздуха в момент времени т в точке с координатами (х, у, ¿) на оси X, У, Z соответственно; Р(х, у, 7,х) - давление воздуха в точке с координатами (х, у, 7); р(х, у, 7,т) - плотность воздуха в точке с координатами (х, у, 7) .

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

/;ггг:с г. ктмжггя I м

-а (77), 201а/ 11 и

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

Для учета турбулентных процессов в пограничном слое используется «подсеточная» полуэмпирическая модель, которая работает только в тех узлах расчетной сетки (или элементах пространства), которые соседствуют с поверхностью обтекаемого тела В настоящей работе рассматривается модель, в которой весь пограничный слой носит турбулентный характер, т е она пригодна для достаточно больших чисел Рейнольдса

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

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

Т = Т, + Т = Т,

Е

фAV

(8)

где Е( - турбулентная энергия воздуха в рассматриваемом объеме А^ с - теплоемкость воздуха; р -плотность воздуха в рассматриваемом объеме

66 №

ЦТ ^ККГТШТГГГГП

4 (77), 2014-

На рис 1 показано реализованное модельное представление пограничного слоя Предполагается, что пограничный слой состоит из движущихся цилиндрических вихревых областей (вихрей), диаметр которых равен толщине слоя 5 . Линейная скорость воздуха, прилегающего к твердой поверхности, равна 0 Воздух на внешнем крае пограничного слоя имеет скорость ламинарного потока их, вычисленную по сеточной модели (Навье-Стокса) . Таким образом, центр вращения вихря (и сам вихрь) перемещается со скоростью их/2 в направлении воздушного потока и вращается с угловой скоростью:

ю = -

и,/2

5/2

Ux 5

(9)

(

т0 = 0,0225pux

1

А

1/4

Res

I 3 7

= 0,02254ПР^Х ,(10)

2тпЛги (ПР3и 11 P = —-- = 0,045Az4 *

(11)

Твердая поверхность Рис . 1 . Модель турбулентного пограничного слоя

Предполагается, что эта механическая энергия затрачивается на поддержание вращения вихря с линейной скоростью наружного края, равной 0./2 . Рассеивание энергии вихря сопровождается ростом его радиуса, что эквивалентно увеличению толщины пограничного слоя

Кинетическая энергия вращающегося вихря:

E =-

J ю2

(12)

где J = mR2 = р52 А7(5 /2)2 - момент инерции вихря; ю - угловая скорость вращения вихря:

E = р52 Azu,2

8

(13)

Согласно экспериментальным исследованиям [6], касательное напряжение турбулентного потока у стенки:

Откуда следует формула для определения толщины пограничного слоя:

5 =

V

8E

PAzUx

(14)

где р - плотность воздуха; ux - скорость набегающего воздушного потока; 5 - толщина пограничного слоя; п - динамическая вязкость воздуха; Re5 = рих5/п - число Рейнольдса .

Механическая мощность этого взаимодействия затрачивается на поддержание вращательного движения вихря:

Pt = M ю ,

T_5Az

где M =-= 2t0Az - вращающий момент, дей-

5 /2 0

ствующий на вихрь со стороны поверхности тела; ю - угловая скорость вращения вихря:

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

cE

dEt ux п

+ P.

dt dx 2

(15)

Уравнение (15) записано для одномерного случая . Его обобщение для двумерного случая:

cE

dEt Ux

dEt Uz

dt dx 2 dz 2

P-

(16)

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

F„ =

X (.P(x,y,z)-P0)AyAz, (17)

—воздух

F3= X (P(x,y,z)-P0)AyAz, (18)

(x,y,z)-B03jiyx (x-Ax,y,z)-TSsio

Z (P(x,y,z)-P0)AxAz, (19)

(л,у,г)-воздух

2

2

2

rvsm/С7

-4 (77), 20ia / U f

X (.P(x,y,z)-P0)AxAz, (20)

(х,у,г)-воздух (*,3>+Áy,z)-TeTO

где Еп, - соответственно сила давления,

действующая на аппарат спереди, сзади, сверху, снизу; Р(х, у, г) - давление в конечно-разностном элементе с координатами (х, у, г); Р0 — 105005 Па -нормальное давление воздуха; Ах, Ау, Аг - шаг конечно-разностной сетки

Силы трения, действующие с разных сторон на аппарат:

7П = X °>0225

(*,;)>,г)-ВОЗДуХ

(х+Дя^^-тено

ЛР 3К(х'У'гУ

5 (x,y,z)

AyAz, (21)

Y3 = X °'0225(

(*,<У,2>-В03ДУХ \

(x-ÁX,y,z)-TOK)

хв = 2 °'0225(

(*,.у,г)-воздух \

(x,y-Áy,z)-Teno

хи = 2 °>0225{

(x,y+Áy,z)-TGm

r]p3Vy(x,y,z)7

5 (x,y,z)

AyAz, (22)

ЛР %(*,y,z)7

8(x,y,z)

AjcAZ, (23)

ripVr^^z)7

SO,}>,z)

AxAz, (24)

= X °'0225

(дг^^-воздух

(25)

б (x,y,z)

Z °'0225^

(*,;)> ,г)-ВОЗДуХ

(*,.)>,z+Az)-iejio

хщ= Z 0>02254

(х,^,г)-воздух (x,_y,z-Az)-icno

г|Р3Г(л:,>>,г)7

AjcAy, (26)

tipXC*'^)7

AxAy, (27)

Хл = Z °>0225

(х,7,г)-воздух (x,_y,z+Az)—ТЕЛО

ЛР 3К(х,У,гУ

8(x,y,z)

AxA_y, (28)

где Уп - сила трения, действующая на аппарат спереди, Уз - сзади, Хв - сверху, Хн - снизу, Упр - справа вверх, Ул - слева вверх, Хпр - справа назад, Хл -слева назад

Лобовое сопротивление:

П 3 в н пр л

Подъемная сила:

Коэффициент лобового сопротивления:

с = 2Х

' ^нРоП2' Коэффициент подъемной силы: 27

Ст =

s„pA2

(29)

(30)

(31)

(32)

где - несущая поверхность; р0 — 1,25 кг/м3 -плотность воздуха на бесконечном удалении от аппарата; У0 — 30 м/с - скорость воздуха на бесконечном удалении от аппарата

Рис . 2 . Распределение давлений вокруг крыльев с профилями серий А12, А15 и А18 при углах атаки набегающего потока 0, 4 и 12°

68 Я

ЦТ ^ККГТШТГГГГП

4 (77), 2014-

Для определения степени адекватности разработанной модели было проведено моделирование обтекания воздухом различных крыльев типовых профилей при разных углах атаки набегающего потока На основе полученных трехмерных полей давлений были рассчитаны основные аэродинамические характеристики крыльев (лобовое сопротивление и подъемная сила), которые были сопоставлены с известными данными [8] по результатам продувки этих крыльев в аэродинамических трубах

На рис 2 показано графическое отображение полей давлений воздуха вокруг крыльев различных профилей при разных углах атаки набегающего потока

На рис 3 приведены вычисленные аэродинамические характеристики соответствующих крыльев в виде поляр, а также поляры, построенные на основании натурных экспериментов

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

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

сопротивления (Сх), коэффициента подъемной силы (Су) и аэродинамическое качество (K).

Процесс моделирования аэродинамических характеристик состоит из следующих этапов:

Э т а п 1. Построение Sd-модели прототипа летательного аппарата (файл планер . stl с помощью SolidWorks) .

Э т а п 2 . Импортирование файла планер . stl в программу «АэроПро-1» и создание сеточной модели

Э т а п 3 . Задание скорости движения летательного аппарата (ио), угла атаки (а), угла скольжения .

Э т а п 4 . Моделирование функций распределения давлений и скорости аэродинамической среды P = fx, у, z, т), и = fx, у, z, т) .

Э т а п 5 . Расчет приложенных к аппарату силовых нагрузок и аэродинамических характеристик Сх, Су, K.

Программа «АэроПро-1» обладает следующими функциональными возможностями:

• расчет распределения скоростей воздуха в аэродинамической среде;

• расчет распределения давлений воздуха в аэродинамической среде;

• получение цифровой информации о значении скорости воздуха в выделенной точке пространства аэродинамической среды для заданного момента времени;

• получение цифровой информации о значении давления воздуха в выделенной точке простран-

Рис . 3 . Сравнение аэродинамических характеристик крыльев серий А12, А15 и А18, полученных модельным и натурным способами

ства аэродинамической среды для заданного момента времени;

• построение графика скорости воздуха в выделенной точке пространства от времени V = ДХ, У, т);

• построение графика давления воздуха в выделенной точке пространства от времени Р = ДХ, У, т);

• визуализация распределения скорости воздуха по сечениям моделируемого пространства;

• визуализация распределения давления воздуха по сечениям моделируемого пространства;

• визуализация распределения скорости воздуха на трехмерной поверхности летательного аппарата;

• визуализация распределения давления воздуха на трехмерной поверхности летательного аппарата;

• трехмерная визуализация скоростных потоков воздуха в аэродинамической среде;

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

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

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

гу.т^г: г: кгткштгта / со

-а (77),201а/ ив

ская среда которого исследована для двух различных углов атаки и заданной скорости В настоящее время разработан математический аппарат для вычисления дополнительных аэродинамических характеристик, необходимых для конструктора планеров, таких, как моменты крена, тангажа и рыскания (при различных отклонениях руля направления, высоты, элеронов и закрылков) На основе этих данных проведено моделирование и расчет группы аэродинамических характеристик, которые позволяют оценить аэродинамические характеристики разрабатываемого планера на этапе проектирования . Развитие этой работы представляет интерес для конструкторов, занимающихся разработкой беспилотных авиационных комплексов Следует отметить, что программа «АэроПро-1» создается на основе собственных математических моделей, которые могут динамично изменяться в зависимости от желаний пользователей

Разрабатываемый инструмент для проектирования планеров позволит существенно снизить время проектирования и помочь конструктору избежать возможных ошибок при проектировании «АэроПро-1» может быть использовано для оценки различных вариантов геометрии планера при предпроектной подготовке и выборе лучших технических решений при создании беспилотных летательных комплексов

Литература

1. Ч и ч к о, А . Н . Клеточно-автоматное моделирование процесса течения расплава в форме / А . Н . Чичко, С . Г Лихоузов // Докл . НАН Беларуси. 2001. Т. 45 . № 4 . С . 110-114 .

2 . Ч и ч к о, А . Н. Комплекс программных средств «ПРОЛИТ» для моделирования процессов течения и охлаждения расплавов / А . Н . Чичко, В . Ф . Соболев, С . Г Лихоузов // Программные продукты и системы . 2002. № 4 .С . 47-48 .

3 . Ч и ч к о, А . Н . Компьютерные системы моделирования физических процессов / А . Н . Чичко // Вестник . 2003 . № 2 . С 42-48

4 . Ч и ч к о, А . Н. Программное обеспечение под суперкомпьютер «СКИФ» «ПроЛит-1с» и «ПроНРС-1с» для литейного и металлургического производства / А . Н . Чичко, Д . М . Кукуй, В . Ф . Соболев, С . Г Лихоузов, Ю . В . Яцкевич, О . И . Чичко // Литье и металлургия. 2008 . № 3 . С . 131-139 .

5 .Г а р б а р у к, А .В . Моделирование турбулентности в расчетах сложных течений: учеб . пособ . / А . В . Гарбарук, М. Х . Стрелец, М . Л . Шур . СПб: Изд-во Политехн . ун-та, 2012 .

6 .Л о й ц я н с к и й Л. Г. Механика жидкости и газа, М. , 6-е изд. 1987.

7 . Ш л и х т и н г Г Теория пограничного слоя / Пер .с нем. М. , 1974 .

8 .К р а в е ц А .С . Характеристики авиационных профилей . М. : Оборонгиз, 1939.