Preview

Вестник Концерна ВКО «Алмаз – Антей»

Расширенный поиск

Исследование влияния параметров метода вихревых элементов на расчет аэродинамических характеристик недеформируемых тел вращения

Полный текст:

Аннотация

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

Для цитирования:


Кудров М.А., Щеглов А.С., Бугаев В.С. Исследование влияния параметров метода вихревых элементов на расчет аэродинамических характеристик недеформируемых тел вращения. Вестник Концерна ВКО «Алмаз – Антей». 2019;(1):51-58.

For citation:


Kudrov M.A., Shcheglov A.S., Bugaev V.S. Study of vortex element method parameters and their effect on rigid rotation bodies aerodynamic computations. Journal of «Almaz – Antey» Air and Space Defence Corporation. 2019;(1):51-58. (In Russ.)

Введение

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

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

Модель вихревого отрезка с конечным радиусом дискретности Вихревой отрезок (рис. 1) представляет собой часть бесконечной вихревой нити и задается следующими величинами:

  • радиус-вектор центра вихревого отрез­ка r0;
  • радиус-вектор конца вихревого отрез­ка r1;
  • интенсивность вихревого отрезка Γ;
  • радиус вихревого отрезка ε.

 

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

 

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

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

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

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

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

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

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

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

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

Вихревые рамки

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

 

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

 

 

Рис. 3. Пример восьмиугольной рамки

 

Четырехугольная рамка задается коорди­натами своих четырех вершин: r01, r11, r21, r31.

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

 

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

 

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

Алгоритм подготовки к расчету

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

  1. Выполнение алгоритма «разбиение тела», после чего имеются совокупность вих­ревых рамок, аппроксимирующих тело, сово­купность контрольных точек Ki, точек вычис­ления давления Kip, нормалей поверхности ni- и площадей рамок Si. Разбиение тела одно­значно определяется числами NL - количество разбиений вдоль образующей, Νφ - количе­ство разбиений по углу относительно оси сим­метрии, np - высота, на которую приподняты точки вычисления давлений относительно кон­трольных точек, а также формой образующей:

Kip = Ki + npni.                                                  (6)

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

  1. Вычисление матрицы А. Элемент ма­трицы Aij = (Qj (Ki), ni), где Q(Ki)- проек­ция вектора скорости, создаваемой j-й рамкой с единичной циркуляцией в контрольной точке Ki на нормаль ni. к этой рамке. При вычисле­нии матрицы для рамок полагается ε = 0.
  1. Вычисление матрицы A-1 один раз в начале расчета.

Алгоритм одного шага расчета

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

  1. Вычисление столбца правых частей системы ΑΓ = b:

bi = (V (Ki), ni),                                                  (7)

где V (Ki) - скорость, индуцируемая всеми свободными вихревыми элементами и набе­гающим потоком в Ki.

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

Γ = A -1b.                                                           (8)

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

П1 - подъем рамок и их расформирова­ние. Образовавшиеся на текущем шаги рамки поднимаются на высоту δup над поверхностью, после чего «разбираются» на составляющие их отрезки и пополняют вихревой след;

П2 - объединение соседних отрезков из рамок. Аналогична процедуре П8 (приведена ниже), но выполняется только для вихревых элементов, рожденных на текущем шаге рас­чета;

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

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

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

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

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

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

к - номер вихревого элемента.

  1. Выполнение следующих процедур:

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

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

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

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

П8 - объединение близких вихревых элементов. Объединяются элементы, расстояние между центрами которых меньше, чем ε* и | cos(a)|> ε**, где α - угол между направляющими векторами h вихревых элементов.

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

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

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

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

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

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

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

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

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

Значение Lmax, при котором влияние вих­ревых элементов на тело можно не учитывать, выбирается, как правило, равным 5-10 харак­терным размерам тела.

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

Результаты расчета

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

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

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

 

Рис. 5. График среднего распределения давления

 

 

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

 

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

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

 

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

 

Исследование влияния параметров про­изводилось следующим образом. Относитель­но вышеуказанных значений изменялось ε в пределах от 0,02 до 0,1 при неизменных зна­чениях остальных величин. На рис. 8 приве­ден график зависимости среднего значения коэффициента сопротивления за последние 200 шагов расчета . Аналогичным образом исследовалось влияние па­раметра nр, который варьировался в пределах от 0 до 0,1 м. На рис. 9 показан график зависимости . Также приведено распределение давления при разной высоте точек вычисления давления над поверхностью пр (рис. 10). На рис. 11 представлено распределение давления при различных значениях параметра ε. На том же графике изображены экспериментальные распределения давления, соответствующие Re = 104 (штриховая линия) и Re = 106 (сплошная линия).

 

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

 

 

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

 

 

 

Заключение

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

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

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

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

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

Об авторах

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


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


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


Для цитирования:


Кудров М.А., Щеглов А.С., Бугаев В.С. Исследование влияния параметров метода вихревых элементов на расчет аэродинамических характеристик недеформируемых тел вращения. Вестник Концерна ВКО «Алмаз – Антей». 2019;(1):51-58.

For citation:


Kudrov M.A., Shcheglov A.S., Bugaev V.S. Study of vortex element method parameters and their effect on rigid rotation bodies aerodynamic computations. Journal of «Almaz – Antey» Air and Space Defence Corporation. 2019;(1):51-58. (In Russ.)

Просмотров: 38


Creative Commons License
Контент доступен под лицензией Creative Commons Attribution 4.0 License.


ISSN 2542-0542 (Print)