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

УДК 533.6:51.001.57

Г. А. Щеглов

МОДИФИКАЦИЯ МЕТОДА ВИХРЕВЫХ ЭЛЕМЕНТОВ ДЛЯ РАСЧЕТА ГИДРОДИНАМИЧЕСКИХ ХАРАКТЕРИСТИК ГЛАДКИХ ТЕЛ

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

Ключевые слова: вихревой метод, вортон, пространственное течение,

гидродинамические характеристики.

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

Для расчета пространственных течений в настоящее время широко применяется метод дискретных вихрей с замкнутыми вихревыми рамками [1]. Он позволяет моделировать завихренность, сосредоточенную в тонких вихревых пеленах, которые сходят в поток по линиям отрыва. Такой подход эффективен при расчете несущих поверхностей и тел с острыми кромками. Однако при расчете обтекания гладких тел возникают значительные трудности при определении линий отрыва. Кроме того, пелена, составленная из связанных между собой вихревых рамок, не допускает "перезамыкания" вихревых линий, что может стать причиной вычислительной неустойчивости.

При расчете обтекания гладких тел используют различные модификации метода вихревых элементов [3, 4], в которых область, занятая завихренностью, моделируется вихревыми элементами - вортонами, которые не связаны между собой и могут моделировать "перезамыкание" вихревых линий. В методе вихревых элементов завихренность генерируется на всей поверхности обтекаемого тела, а области отрыва потока формируются при расчете естественным образом. Недостатком метода вихревых элементов является нарушение соленоидально-сти поля завихренности при аппроксимации его вихревыми элементами [5]. Поэтому при расчетах необходимо использовать большое число вортонов и контролировать соленоидальность, что требует дополнительных вычислительных ресурсов [3]. Для повышения эффективности метода проводится поиск алгоритмов аппроксимации поля завихренности малым числом вихревых элементов с достаточной для инженерной практики точностью. Также проводится поиск наилучшей расчетной схемы для моделирования процесса генерации завихренности вблизи поверхности тела.

В предыдущих работах автором было проведено успешное тестирование вихревым методом симметричного вортона-отрезка [6] и рамок, составленных из вортонов, для формирования расчетной схемы на поверхности тела [7].

Цель настоящей работы — объединение указанных моделей для расчета гидродинамических характеристик гладких тел и тестирование полученной модификации метода вихревых элементов.

Постановка задачи и уравнения движения. Исследуется пространственное обтекание тела мгновенно стартующим потоком несжимаемой среды плотностью рз = const при числах Рейнольдса порядка 105. В этом случае возможно применение подхода Прандтля. Влияние вязкости учитывается как причина генерации завихренности в тонком пристеночном слое вблизи поверхности тела, а движение среды во всем остальном пространстве рассматривается в рамках модели идеальной жидкости. Обоснование такого подхода применительно к движению завихренного потока дано в работе [1]. Движение среды описывают уравнениями неразрывности и Гельмгольца:

div V = 0;

dt = rot(V X й), (1)

где й (r, t) - нестационарное трехмерное поле скоростей; r - радиус-вектор точки в неподвижной мировой системе координат; й = rotV -завихренность.

Задано граничное условие затухания возмущений на бесконечности

lim й (r, t) = йзо = const

и условие непротекания на поверхности тела £:

V ,г) х из = 0, (2)

где гз - радиус-вектор точки на поверхности тела, заданный в мировой системе координат;

- единичный вектор внешней нормали к поверхности тела в этой точке.

При дискретизации уравнений движения жидкости завихренность аппроксимируется множеством из N вихревых элементов (ВЭ) - симметричных вортонов-отрезков (рис. 1). Такой ВЭ подобен элементу вихревой трубки. Он показал высокую эффективность при моделировании явления "чехарды" вихревых колец [6].

Эволюция поля завихренности складывается из движения центров ВЭ по жидким линиям, изменения длины, поворота направляющих векторов ВЭ (при движении интенсивность ВЭ не изменяется) и генерации новых ВЭ.

Уравнения движения для центра и направляющего вектора г-го ВЭ имеют вид

Рис. 1. Симметричный вортон-отрезок:

г0 — радиус-вектор центра в мировой системе координат; Дг — направляющий вектор

drr

0i

dt

= V (roi, t);

d (Ari) dt

N

(r0i) ) x Ari, i = 1,..., N,

vj = 1

(3)

(4)

где

N

V (r0i, t) = + ^ Vj (roi); Vj (rOi) = (4п)-1 Гc a; j=i

c = 2 ((771 - 72) x Arj)/h; a = 2 (So x Arj); h = |a|2; So = r0j - rOi; 77I = si/si; П2 = S2/S2; = So + Arj; S2 = So - Arj; si = |Si|; S2 = |S21; Ar2 = |Arj|2;

Bj штп =

dc

dc

d (7 )m = Г d (7o)n 4n Vd (7oi) 2

x (a)m - 2 cemnk (A7j)

j>k / '

/ dh (h (Si - S2) (A7j)n - дрТ^K) +

d (7o)n 2

+

h2Si S2

h SiS2

d (7oi)i

((A7j x Si) s2 ni - (A7j x S2) Si 772) ;

n

dh

:-8 • ((So)n Ar? - (Arj)n X (Arj • So)) ;

d (rbi )n v

K = (si + S2) Arf - (si - S2) • (Afj x 50) ;

m,n,k = x, y, z - компоненты векторов; £mnk — тензор Леви-Чивиты.

Для исключения неограниченного роста скоростей Vj (foi ) и их производных при приближении к оси ВЭ, вводится радиус вихревой трубки е, т.е. при |a| / (2 |Afj |) < £ индуцированные скорости и их производные убывают по линейному закону до нуля на оси ВЭ:

V; (foi) = aVj (fOi) ; в; (f*)mn = a Bj (f0i)mn ,

где fi, = foi + 5 (a-i - 1), a = |3|/ (2е |Af;|), f = Af;So X Af - *o.

A j

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

Генерация завихренности моделируется рождением ВЭ из M вор-тонных рамок, образующих расчетную схему на поверхности тела (рис. 2). Вортонные рамки обеспечивают замкнутость и соленоидаль-ность вихревой системы вблизи поверхности. Такая схема была рассмотрена ранее с использованием сферических вортонов [7]. Каждая j-я m-угольная вихревая рамка задана радиусами-векторами вершин fS (s = 1,...,m), циркуляцией Г;, радиусом-вектором контрольной точки k0 и внешней нормалью П0 к поверхности в контрольной точке.

Вихревые отрезки, составляющие рамку, являются симметричными вортонами-отрезками с параметрами:

; = 1 (fS + fj+i ) + e no ; Aff = 1 (fj+i - fj ) ; (5)

rf = rj; s = 1,..., m, 0 < в < е.

Рис. 2. Вортонные рамки на поверхности обтекаемого тела

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

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

М т

Е Г1 Е X П0 = X П0, з = 1,..., М, (6)

¿=1 ^=1

где г7в (г) = К (г) — вектор влияния вортона с единичной интен-

Г3 = 1

сивностью. Корректное решение уравнения (6) должно обеспечивать нулевую скорость в точках, лежащих внутри тела. Поскольку, как показано в работе [1], для замкнутой поверхности система вырождена (6), вводится дополнительное условие

М

=0, з = 1,...,м,

1

и регуляризирующая переменная. После решения уравнения (6) происходит разделение рамки на ВЭ, которые пополняют вихревой след.

Давление в контрольных точках рамок определяется при помощи аналога интеграла Коши-Лагранжа [8]:

p (r, t) = p+ p^

у2 у (Г* t)2 ^

- + ЕУ (ri,t) X у (r) - Ib (t)

2 2

i=1

(7)

где — интеграл, учитывающий генерацию вортонов; p^ = const — давление в невозмущенном потоке. Главный вектор сил, действующий на обтекаемое тело, находится как

м

R (t) = - ЕpfU^jO,

j=1

где Sj - площадь панели j-й рамки. Гидродинамические коэффициенты определяют по известным формулам.

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

Алгоритм решения задачи. Алгоритм расчета одного шага состоит из следующих операций.

1. Построение расчетной схемы из M вортонных рамок на поверхности тела и определение их параметров (5).

2. Решение системы (6) и нахождение неизвестных циркуляций Г (j = 1,...,M ).

3. Пополнение вихревого следа новыми ВЭ. При этом элементы, расположенные на сторонах соседних рамок, объединены в один ВЭ (процедура объединения приведена далее в п. 6). Также для уменьшения числа ВЭ в следе элементы с малой интенсивностью |Щ < rmin (rmin - заданная малая константа) не пополняют след.

4. Определение давления в контрольных точках по формуле (7) и вычисление гидродинамических нагрузок и коэффициентов.

5. Численное интегрирование уравнений эволюции вихревого следа (3), (4) методом первого порядка с шагом At и определение параметров ВЭ на новом шаге. При этом ВЭ, попавшие внутрь тела, возвращаются в поток. Для повышения устойчивости счета проводится контроль приращения векторов An Если приращение AF за один шаг расчета превышает заданную величину ед, то приращение направляющего вектора ВЭ обнуляется.

6. Реструктуризация ВЭ в следе, т.е. объединение ВЭ и ограничение длины вихревого следа. Объединение элементов, для которых

Г • • Г '

выполнены условия: |n — rj| < е* и |Ф| < е**, где Ф = j (если

N |rj1

Ф < 0, то производится разворот вектора вортона Arj с изменением знака циркуляции Г ). Параметры объединенного ВЭ задаются следующими выражениями:

|Г»| r + |Г,-| rj . |Г»| Ari + |Г,-| Arj _ _ _

rs =1 ;Г: + 'jI j ; Ars =1 г| lr: +1 j|—j; rs = r + | г • 1 + |rj 1 | Г •| + |rj 1

Цикл вычислений продолжается до достижения конечного времени вычислений ^ или до остановки вычислений пользователем. При реализации алгоритма разработана программа MVE3D, в которой использована библиотека MPI для распараллеливания 5-й и 6-й операций. Расчеты проводились с использованием кластера МВС-100к.

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

В численных экспериментах дискретизация поверхности тел вращения вортонными рамками характеризовалась двумя параметрами: N - число рамок по образующей тела вращения; N2 - число рамок в окружном направлении. Все рамки являются четырехугольными, за исключением полюсов, где использовались многоугольные рамки с числом сторон N2. Таким образом, число вортонных рамок M = = N1 N2 + 2.

В каждом вычислительном эксперименте время окончания счета равнялось £к = 80, а шаг ДЬ = 0,04, т.е. проводилось 2000 шагов расчета. Нестационарные коэффициенты сопротивления осреднялись по 1500 временных шагов (на интервале от 500-го до 2000-го шага). Полученные таким образом стационарные коэффициенты сопротивления сравнивались с экспериментальными данными [9].

На первом этапе исследований выбирали безразмерные параметры расчетной схемы. В качестве тестового объекта был выбран тонкий диск диаметром Б = 1 и толщиной Н = 0,1. В диапазоне чисел Рей-нольдса Яе от 103 до 107 стационарный коэффициент сопротивления для характерной площади 5 = 0,25 пБ2 имеет постоянное значение СХхр = 1,17 [9].

По результатам вычислительных экспериментов были выбраны безразмерные параметры = {1, 0, 0}т; рж = 1; р^ = 1; в = 0,01; г = 0,1; г* = 0,7 г; г** = 0,9; гд = 0,02; Г^ = 10-6; £/аг = 10,0.

В качестве примера на рис. 3 приведена зависимость изменения коэффициентов СХа, Суа, CZa от времени для диска с параметрами N = 22, И2 = 36. На рис. 4 показан вид вихревого следа, сформированного N = 8667 ВЭ (отмечены на рисунке точками), на последнем шаге расчета (Ь = 80).

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

На втором этапе исследований, используя определенные ранее безразмерные параметры, проводили моделирование обтекания шара.

Рис.3. Зависимость коэффициентов Сх а, Су а, Сг а от времени для диска (штриховая — экспериментальное значение)

б

Рис. 4. Вихревой след за диском (а) и шаром (б)

В расчетах [7] был получен стационарный коэффициент сопротивления шара, близкий к докритическому значению, равному 0,47. С использованием симметричных вортонов-отрезков для шара с параметрами N = 10, N2 = 24 получен коэффициент СХа = 0,101, для шара с параметрами N1 = 17, N = 34 - СХа = 0,252, для шара с параметрами N = 20, N = 40 - СХа = 0,224. Таким образом, с увеличением дискретизации вычисленный коэффициент сопротивления приближается к закритическому значению 0,2. На рис.4,б в качестве примера показан вихревой след, сформированной N = 3763 ВЭ (отмечены на рисунке точками), в конце расчета (£ = 80).

Таблица

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

n х n2 Вычисленный CX Погрешность, %

10 х 24 0,87 25,3

22 х 24 1,05 10,2

22 х 30 1,08 7,3

22 х 36 1,13 3,1

40 х 38 1,18 1,8

Поскольку кризис обтекания шара происходит в диапазоне чисел Рейнольдса от 105 до 4 • 105, можно сделать вывод, что для моделирования течений среды с малой вязкостью предпочтительнее использовать симметричный вортон-отрезок, чем ВЭ из работ [3, 4].

Рис. 5. Сравнение расчетного (точки ■ и А для двух расчетных схем с различными 1 х N2) = 24 х 16 и 36 х 30) и экспериментального (сплошная кривая) [9] коэффициентов СХа (а) и Суа (б) для цилиндра

На третьем этапе исследований рассчитывали стационарные коэффициенты СХа и Суа кругового цилиндра диаметром Б = 1 и длиной Ь = 2 в зависимости от угла атаки. Экспериментальные значения стационарных коэффициентов для характерной площади 5 = ЬБ приведены в работе [9].

Расчеты проводили с безразмерными параметрами, определенными на первом этапе исследований. На рис. 5, а и б показаны расчетные (точки для двух расчетных схем с различными N х N) и экспериментальные значения для коэффициентов СХа и Суа.

Как видно из рис. 5, для коэффициента СХа наблюдается хорошее совпадение расчетных и экспериментальных данных. В то же время для коэффициента Суа на углах атаки от 3° до 50о при расчете получаются завышенные значения. Анализ расчетных данных показал, что хотя в этих случаях решение системы (6) оказывается корректным: условие непротекания в контрольных точках выполняется, скорость внутри тела, вычисленная и усредненная по 100 точкам, оказывается отличной от нуля, составляя более 20% от | вне зависимости от дискретизации тела. В то же время контроль скорости внутри тела после решения уравнения (6) для диска и сферы дает значение не более 3 % от |, которое уменьшается с увеличением дискретизации.

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

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

инженерных расчетов точностью. С увеличением числа вортонных рамок на теле точность определения гидродинамических характеристик повышается.

Алгоритм метода вихревых элементов может быть эффективно распараллелен. Время выполнения 2 000 временных шагов при расчете на одном ядре процессора Intel Core2Duo 2,6 ГГц занимал 8 ч, при параллельном расчете на двух ядрах время вычислений снижалось практически вдвое. С увеличением числа использованных процессоров эффективность распараллеливания снижалась. При использовании 32 процессоров кластера МВС-100к тот же расчет занимал около 30 мин.

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

Автор благодарит Межведомственный суперкомпьютерный центр РАН за предоставленную возможность использования высокопроизводительного кластера МВС-100к.

СПИСОК ЛИТЕРАТУРЫ

1. Трехмерное отрывное обтекание тел произвольной формы / Под ред. С.М. Бе-лоцерковского. - М.:ЦАГИ, 2000. - 265 с.

2. L e w i s R. I. Vortex Element Methods for Fluid Dynamic Analysis of Engineering Systems // Cambridge University Press, 2005. - 567 p.

3. PloumhansP., WinckelmansG. S., Salmon J. K., Leonard A., Warren M. S. Vortex Methods for Direct Numerical Simulation of Three-Dimensional Bluff Body Flows: Applocation to the Sphere at Re=300, 500 and 1000 // Journal of comp. Physics. - 178. - 2002. - P. 427-463.

4. O j i m a A., K a m e m o t o K. Numerical simulation of unsteady flow around three dimensional bluff bodies by an advanced vortex method // JSME International Journal, Series B. - Vol. 43. - No. 2 (2000). - P. 127-135.

5. К о р н е в Н. В. Метод вихревых частиц и его приложение к задачам гидроаэродинамики корабля / Дис.... д-ра техн. наук. - СПб, 1998. - 184 с.

6. Марчевский И. К., Щеглов Г. А. Модель симметричного вортона-отрезка для численного моделирования пространственных течений идеальной несжимаемой среды // Вестник МГТУ им. Н.Э. Баумана. Серия "Естественные науки". - 2008. - № 4. - С. 62-71.

7. Щеглов Г. А. О применении вортонных рамок в методе вихревых частиц // Вестник МГТУ им. Н.Э. Баумана. Серия "Естественные науки". - 2008. - № 2. - С.15-22.

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

9. Д е в н и н С. И. Аэрогидромеханика плохообтекаемых конструкций. - Л.: Судостроение, 1983. - 320 с.

Статья поступила в редакцию 23.06.2008

Георгий Александрович Щеглов родился в 1972 г., окончил МГТУ им. Н.Э. Баумана в 1996 г.. Канд. физ.-мат. наук, доцент кафедры "Аэрокосмические системы" МГТУ им. Н.Э. Баумана. Автор более 30 научных работ в области динамики и прочности конструкций аэрокосмических систем, вычислительных методов в гидроупругости.

G.A. Shcheglov (b. 1972) graduated from the Bauman Moscow State Technical University in 1996. Ph. D. (Phys.-Math.), assoc. professor of "Aerospace Systems" department of the Bauman Moscow State Technical University. Author of more than 30 publications in the field of dynamics and strength of structures of aerospace systems, computation methods in hydro-elasticity.