Исследование влияния параметров метода вихревых элементов на расчет аэродинамических характеристик недеформируемых тел вращения Текст научной статьи по специальности «Физика»

К<1

УДК 532.527

М. А. Кудров, А. С. Щеглов, В. С. Бугаев Исследование влияния параметров метода вихревых элементов на расчет аэродинамических характеристик недеформируемых тел вращения

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

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

Введение

Для расчета нестационарных аэродинамических характеристик тел при малых числах Маха (М <0,3.„0,4), когда сжимаемостью среды можно пренебречь, целесообразно использование бессеточных методов с применением вихревых элементов. Практическая ценность таких методов состоит в высокой скорости проведения расчетов относительно сеточных методов. Дополнительным преимуществом является отсутствие расчетной сетки, построение которой - непростая задача при моделировании движения тел.

Цель использования данного метода -расчет аэродинамических характеристик тел при больших числах Re ~ 105...106, так как вязкость в нем считается малой и учитывается лишь как причина образования завихренности на поверхности обтекаемого тела. В остальной области пространства применяют модель идеальной жидкости. Научной основой данного метода является теорема Томсона для идеальной жидкости, которая гласит, что циркуляция скорости по произвольному жидкому контуру сохраняется во времени. Это означает, что в области течения завихренность не образуется. Ее появление возможно лишь на поверхности обтекаемого тела. Далее эта завихренность становится свободной. Модель вихревого отрезка с конечным радиусом дискретности Вихревой отрезок (рис. 1) представляет собой часть бесконечной вихревой нити и задается следующими величинами:

© Кудров М. А., Щеглов А. С., Бугаев В. С., 2019

ка Г

0'

Рис. 1. Модель вихревого отрезка радиус-вектор центра вихревого отрез-

радиус-вектор конца вихревого отрез-

ка г,

• интенсивность вихревого отрезка Г;

• радиус вихревого отрезка е.

Поле скорости, индуцируемое вихревым отрезком с нулевым радиусом в произвольной точке пространства г, может быть найдено по формуле Био - Савара [1]:

V = — ca, 4п

(1)

где c = a

1-2

,s2

VI 2I

h = Г1 - Г0

s0 = Г - г0

^ = s 0 - Ъ

S2 = So + ^

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

е

о р

ст о

оте

к

а р

а

m

о Ч е л с с

к с е

у

и м с о К

V (г) =

Г

—су а у, Л >г 4п

Г Я

--сЕаЕ, Я <г,

4п г

0, Я = 0,

N

где а Е= h (гЕ- го);

с. = а .

^ = г -

1-2

_ гв~ го - h Л

|гв" г0 + h| к_ г0 _ h|

^-1 R

в0 h

ь2

Я = 1

(2)

£ = V(г)■

М

(3)

Аналогично В(г0) = £ В. (г0), где В. (г0) -

.=1

11

Уравнения движения вихревого отрезка

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

тензор деформации I - го вихревого отрезка. Для компонент тензора В. (г) от одного вихревого отрезка в произвольной точке пространства имеется аналитическая формула:

В™ (г ) =

Г

4п

где

дРп

= г

^) + ^ Я

дРпУу)т у дрп ) Я

-втп (г,), Я <е,

0, Я = 0,

ткп )k ;

дсу дРп

■ = а„

О - О + ^Су (^ h ))(Ь )п -

V 182| |81| у1

^п »2 (^п «1

V I 21

"2су (80 )п h

(5)

о см

<

I

(0 га

г

о ^

со га г о.

<и

о

и <и со

см ■ч-ю о

I

см ■ч-ю см

(П (П

Записав это уравнение применительно к центру отрезка г0 и к концу отрезка г1, можно получить систему уравнений, описывающую динамику движения вихревого отрезка с точностью о(Ь):

в 00) =

М- = V(Го, *)

м

МЬ о/ чи — = в(г0)h М

дУ дУ дУ

у(г0) уЫ у(г0)

дх ду дz

(4)

дУу _у

дх

дУу_ ду

^ (г0) ^ Ос) ^ (г0)

дУу_

дz

дУ,

дУ^

У

(г0) —^ (г0) —^ (г0)

дх

ду

д

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

Вихревая рамка представляет собой замкнутую систему вихревых отрезков. В работе используются два вида рамок (рис. 2, 3).

Четырехугольная рамка задается координатами своих четырех вершин: г01, г11, г21, г31.

Правильная многоугольная рамка характеризуется координатами своего центра г0, ко -

где В(г0) - тензор деформации вихревого отрезка [2].

V (г) - поле скорости в точке г, которое является векторной суммой скорости потока ^ и полей скорости, индуцируемых всеми N вихревыми элементами, имеющимися в потоке

£ V., т. е. V (гз) = V. (г,).

.=1 .=1

Рис. 2. Пример четырехугольной рамки

в

ординатами двух своих вершин г01, г1Х и числом N отрезков, составляющих рамку. Разбиение тела на рамки Для расчета обтекаемое тело аппроксимируется многогранником, поверхность которого состоит из вихревых рамок двух вышеупомянутых типов. В обоих полюсах тела располагаются правильные многоугольные рамки. Остальные участки покрываются с помощью четырехугольных рамок (рис. 4). Точками обозначены центры вихревых отрезков, сошедших в поток с поверхности тела.

Рис. 4. Пример разбиения сферы с помощью вихревых рамок

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

Подготовка к расчету совершается согласно представленному ниже алгоритму.

1. Выполнение алгоритма «разбиение тела», после чего имеются совокупность вихревых рамок, аппроксимирующих тело, совокупность контрольных точек К., точек вычисления давления Кр, нормалей поверхности п. и площадей рамок Sl. Разбиение тела однозначно определяется числами NL - количество разбиений вдоль образующей, Nv - количество разбиений по углу относительно оси симметрии, пр - высота, на которую приподняты точки вычисления давлений относительно контрольных точек, а также формой образующей:

Кр = К. + прп. (6)

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

2. Вычисление матрицы А. Элемент матрицы А = (Qj (К), п), где Q} (К) - проекция вектора скорости, создаваемой у-й рамкой с единичной циркуляцией в контрольной точке К на нормаль п;. к этой рамке. При вычислении матрицы для рамок полагается е = 0.

3. Вычисление матрицы А-1 один раз в начале расчета.

Алгоритм одного шага расчета Один шага расчета выполняется в представленной ниже последовательности.

1. Вычисление столбца правых частей системы АГ = Ь:

ь = (V (К {),п), (7)

где V (К(.) - скорость, индуцируемая всеми

свободными вихревыми элементами и набе- ~

гающим потоком в К. |

2. Нахождение циркуляций рамок, & рожденных на данном шаге на поверхности о тела, по формуле: £

Г = А -1Ь. (8) *

I

3. Выполнение процедур: 8

П1 - подъем рамок и их расформирова- ¡5

ние. Образовавшиеся на текущем шаги рамки о

поднимаются на высоту Ъир над поверхностью, *

после чего «разбираются» на составляющие их |

отрезки и пополняют вихревой след; ^

П2 - объединение соседних отрезков из у рамок. Аналогична процедуре П8 (приведена

й!

о см

<1

I

м га

г |

0 ^

со га

1

о.

<и

о

и <и со

см ■ч-ю о

I

см ■ч-ю см

(П (П

ниже), но выполняется только для вихревых элементов, рожденных на текущем шаге расчета;

П3 - уничтожение отрезков малой интенсивности. Если существует вихревой элемент со значением интенсивности | Г |< Гт1п, то такой вихревой элемент удаляется из расчета.

4. Расчет давлений в контрольных точках и вычисление гидродинамических сил и коэффициентов.

Давление в произвольной точке пространства г вычисляется с помощью обобщенного интеграла Коши - Лагранжа [3]. Сила, действующая на тело, рассчитывается путем суммирования сил, действующих на каждую панель тела:

R = -£ р(Кр )Si п. ■ (9)

I

5. Численное интегрирование системы (4) методом Эйлера первого порядка с шагом т.

Для каждого вихревого элемента вычисляется новое положение его центра г0 и новое значение направляющего вектора Ь:

ОХ1 = (г,)п + V (г0)к т;

Ь п+1 = Ь к + В(г0)к Ь ? т,

(10)

щего вектора Ь относительно оси за один шаг превышает фтах, то приращение направляющего вектора на данном шаге обнуляется;

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

П8 - объединение близких вихревых

элементов. Объединяются элементы, расстоя-

*

ние между центрами которых меньше, чем е и | cos(a)|> е**, где а - угол между направляющими векторами Ь вихревых элементов.

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

ШК,),п,) ШК,),п,) ... ^(К,), п,) 1 (Ql(K 2), П 2) (Q2(K2),n2) ... ^(К2), П2) 1

... ... ... ... 1

(Ql(K N ), П N ) ), П N ) ... ^ (К N ), nN ) 1

1 1110

(11)

где п - номер шага расчета;

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

П5 - удаление вихревых элементов, расположенных далеко от тела. Если расстояние от некоторого вихревого элемента до центра тела превышает Цагх, то вихревой элемент удаляется из расчета, так как его влияние на тело незначительно;

П6 - контроль приращения Ь и поворота вихревого элемента. Если для некоторого вихревого элемента приращение длины направляющего вектора Ь за один шаг превышает ед, то приращение направляющего вектора на данном шаге обнуляется. Если для некоторого вихревого элемента угол поворота направляю-

(Г1 > ( (V(К1),П1) )

Г2 (V (К 2 ),П2 )

ГN (V(КN)

1 Я ) 1 0 J

Таким образом, результаты расчета определяются следующими параметрами:

NL, Nф, е> Т еА' Фтах, е*> е", п", 5 up, Н, Lmаx, Гтт ■

Методика выбора параметров расчета

Исходя из проведенных вычислительных экспериментов и их сопоставления с экспериментальными данными, была составлена методика подбора параметров.

NL является свободным параметром, а Nф подбирается таким образом, чтобы размеры получающихся при разбиении тела рамок были как можно более одинаковы.

Параметр е выбирается, как правило, равным примерно половине среднего разме-

X

X

l<<«;

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

Шаг интегрирования т должен быть сопоставим со временем прохождения вихревым элементом расстояния, равного характерному

т L

размеру рамки L, т. е. т .

ж

Поворот вихревого отрезка и приращение его длины за один шаг интегрирования

не должен быть слишком большим. Значения

п —

9max = — и ед=0,05L показали хорошие результаты в модельной задаче перезамыкания вихревых колец.

Значение параметра разворота H составляет примерно половину характерного размера рамки. Разворот вихревых отрезков вблизи тела предотвращает возникновение циркуляции поля скорости внутри тела.

Значение Lmax, при котором влияние вихревых элементов на тело можно не учитывать, выбирается, как правило, равным 5-10 характерным размерам тела. Разработка программного комплекса Программный комплекс, реализующий данный алгоритм, выполнен с использованием кроссплатформенного фреймворка Qt, что обеспечивает его функционирование в большинстве существующих операционных систем. В созданном программном комплексе параллельные вычисления осуществляется с помощью высокоуровневого API QtConcurrent. Картина расчета визуализируется с помощью библиотеки OpenGL в режиме реального времени и представляет собой 3D-отображение сетки из рамок на поверхности тела и вихревых элементов. Результаты расчета сохраняются в формате, пригодном для обработки постпроцессором ParaView. Результаты расчета

Расчет останавливается по достижении заданного числа итераций. Число итераций было принято равным Nsteps = 400 для всех расчетов.

В ходе работы исследовалась зависимость значений коэффициента сопротивле-

ния сферы Сх, а также распределения давления вдоль образующей Ср (9) в зависимости от основных параметров расчета. Сравнение расчетных данных с экспериментальными осуществлялось с использованием книги [4], в которой приведены экспериментальные данные стационарных значений коэффициента Сх сферы для различных чисел Re.

На рис. 5 приведен график среднего распределения давления Ср (9) за последние 200 шагов, а на рис. 6 отображены значения Сх за тот же временной промежуток. Расчет был проведен со следующими значениями параметров:

= 50, = 40, е = 0,04 м, т = 0,0785 с, ед = 0,003, фтах = Ц, а* = 0,028 м, а** = 0,9, пр = 0,039 м, 8нр = 0,00785 м, Н = 0,0785 м,

Lmax = Юм, Гт1п = 10-6 м2/с.

Средняя длина панели для такого разбиения оказалась равной АЬ = 0,06 м .

"X

\

\ у

\

\

\ У

> /

/

0,8 0,6 0,4 0,2 0

-0,2 -0,4 -0,6

0 0,5 1,0 1,5 2,0 2,5 0,рад Рис. 5. График среднего распределения давления Сг

0,40 0,35 0,30 0,25 0,20 0,15 0,10 0,05 0

200

240

280

320

360

N

Рис. 6. Коэффициент сопротивления в зависимости от номера шага расчета

е

о р

ст о

оте

к

а р

а

m

о ч е л с с

к с е

у

и м с о К

Среднее значение коэффициента сопротивления за последние 200 шагов оказалось равным Сх = 0,33. Из сопоставления с результатами эксперимента [4] можно сделать вывод, что обтекание с данными параметрами соответствует числу Рейнольдса Re = (2...5)-105. На рис. 7 приведена визуализация картины вихревого следа на последнем шаге расчета. К этому моменту расчета в следе было 18 876 вихревых элементов.

0,4 0,3 0,2 ОД

Сх(пр

)

0 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6 пр

Рис. 9. Среднее значение коэффициента сопротивления от высоты подъема точек вычисления давлений

о см

<

I

о га

г |

о ^

со га г о.

<и

о

и <и со

см ■ч-ю о

I

см ■ч-ю см

(П (П

Рис. 7. Картина вихревого следа на последнем шаге расчета

Исследование влияния параметров производилось следующим образом. Относительно вышеуказанных значений изменялось е в пределах от 0,02 до 0,1 при неизменных значениях остальных величин. На рис. 8 приведен график зависимости среднего значения

коэффициента сопротивления за последние

— _ _ ^

200 шагов расчета Сх(е), где е=—■ Анало-

Т

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

от 0 до 0,1 м. На рис. 9 показан график зави-

- _ _ пр

симости Сх(пр), где пр= —. Также приведе-

А Т

но распределение давления при разной высоте

точек вычисления давления над поверхностью пр (рис. 10). На рис. 11 представлено распре-

0,6 0,5 0,4 0,3 0,2 ОД

СЛ Ё)

0 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6 г

Рис. 8. Среднее значение коэффициента сопротивления в зависимости от е

0,8 0,6 0,4 0,2 0 -0,2 -0,4 -0,6

хоА

г А

«Д г*

ощ 0> (0) ,0 ллллллл

с йуаш ш

о "*х> О °о<х х*§8 0°

0

0,5 1,0 1,5 2,0 2,5

0, рад

Рис. 10. Распределение давления для различных значений параметра пр:

Л - пр = 0,02 м; О - пр = 0,04 м; * - пр = 0,1 м

0,5 0

-0,5 -1,0 -1,5

'5х 1?о*хххх> Л^ооо« Xх „о»о ?Ч8рроо уКххххх «хххх** ■ **

\ с ,(0) /

0 0,5 1,0 1,5 2,0 2,5 3,0 9, рад

Рис. 11. Распределение давления для различных значений параметра е : Л - £ = 0,04 м; О - е = 0,06 м; « - £ = 0,08 м;

___- Re = 10 (экспериментальные данные);

— - Re = 106 (экспериментальные данные)

l<<«;

деление давления при различных значениях параметра е. На том же графике изображены экспериментальные распределения давления, соответствующие Re = 104 (штриховая линия) и Re = 106 (сплошная линия). Заключение

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

Достоинства метода вихревых элементов:

• гораздо более высокая производительность по сравнению с сеточными методами;

• возможность распараллеливания большинства процедур метода;

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

• отсутствие схемной вязкости.

Недостатки метода вихревых элементов:

• ограниченная область применения по числу Маха ( M <0,3...0,4 );

• несколько более низкая точность результатов расчетов по сравнению с сеточными методами.

Список литературы

1. Щеглов Г. А. Модификация метода вихревых элементов для расчета аэродинамических характеристик гладких тел // Вестник МГТУ им. Н. Э. Баумана. Сер. Машиностроение. 2009. № 2. С. 26-36.

2. Marchevsky I. K., Scheglov G. A. 3D Vortex structures dynamics simulation using vortex fragmentons. URL: https://www.researchgate. net/publication/28693 95 52_3D_vortex_ structures_dynamics_simulation_using_vortex_ fragmentons (дата обращения 10.12.2018).

3. Андронов П. Р., Гувернюк С. В., Дыннико-ва Г. Я. Вихревые методы расчета нестационарных гидродинамических нагрузок. М.: Изд-во МГУ, 2006. 184 с.

4. Виноградов Н. И, Крейндель С. А, Лев И. Г., Нисневич М. З. Привязные подводные системы: Аэродинамические характеристики при установившемся движении. СПб.: Изд-во СПбГУ, 2005. 304 с.

Поступила 01.03.19

Кудров Максим Александрович - кандидат технических наук Федерального государственного автономного образовательного учреждения высшего образования «Московский физико-технический институт (государственный университет)» факультета аэромеханики и летательной техники, г. Жуковский.

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

Щеглов Андрей Сергеевич - инженер Федерального государственного автономного образовательного учреждения высшего образования «Московский физико-технический институт (государственный университет)» факультета —

О

аэромеханики и летательной техники, г. Жуковский. |

Область научных интересов: динамика полета, вихревые методы расчета аэродинамических характеристик лета- о

о.

тельных аппаратов. ^

0

1

Бугаев Василий Сергеевич - техник Федерального государственного автономного образовательного учреждения те.

высшего образования «Московский физико-технический институт (государственный университет)» факультета з

аэромеханики и летательной техники, г. Жуковский. 1

Область научных интересов: область научных интересов: численное решение прикладных задач аэродинамики и §

О

динамики полета летательных аппаратов. ц

Study of vortex element method parameters and their effect on rigid rotation bodies aerodynamic computations

The purpose of this study was to implement a software package that allows non-stationary aerodynamic computations of fixed rotation bodies using the vortex element method. In the course of the work, we developed an algorithm for rigid rotation bodies aerodynamic computations by means of this method. Furthermore, we studied the influence of calculation parameters on the results obtained.

Keywords: vortex segment, vortex frame, Cauchy - Lagrange integral, drag coefficient, pressure distribution, deformation tensor.

Kudrov Maksim Aleksandrovich - Candidate of Engineering Sciences, Federal State Autonomous Educational Institution of Higher Education "Moscow Institute of Physics and Technology (State University)", Department of Aeromechanics and Flight Engineering, Zhukovsky.

Science research interests: development of software system for computer simulation of physical processes.

Shcheglov Andrei Sergeevich - engineer, Federal State Autonomous Educational Institution of Higher Education "Moscow Institute of Physics and Technology (State University)", Department of Aeromechanics and Flight Engineering, Zhukovsky. Science research interests: flight dynamics, vortex methods for aerodynamic computations.

Bugaev Vasiliy Sergeevich - technician, Federal State Autonomous Educational Institution of Higher Education "Moscow Institute of Physics and Technology (State University)", Department of Aeromechanics and Flight Engineering, Zhukovsky. Science research interests: numerical solution of applied aerodynamics and flight dynamics problems.