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

Введение

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Четырехугольная рамка задается коорди­натами своих четырех вершин: r₀1, r₁₁, r₂₁, r₃₁.

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

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

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

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

1. Выполнение алгоритма «разбиение тела», после чего имеются совокупность вих­ревых рамок, аппроксимирующих тело, сово­купность контрольных точек K_i, точек вычис­ления давления K_i^p, нормалей поверхности n_i- и площадей рамок S_i. Разбиение тела одно­значно определяется числами N_L - количество разбиений вдоль образующей, Ν_φ - количе­ство разбиений по углу относительно оси сим­метрии, n^p - высота, на которую приподняты точки вычисления давлений относительно кон­трольных точек, а также формой образующей:

K_i^p = K_i + n^pn_i.                                                  (6)

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

2. Вычисление матрицы А. Элемент ма­трицы A_ij = (Q_j (K_i), n_i), где Q_j (K_i)- проек­ция вектора скорости, создаваемой j-й рамкой с единичной циркуляцией в контрольной точке K_i на нормаль n_i. к этой рамке. При вычисле­нии матрицы для рамок полагается ε = 0.

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

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

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

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

b_i = (V (K_i), n_i),                                                  (7)

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

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

Γ = A ^-1b.                                                           (8)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Заключение

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

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

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

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

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