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

УДК 533.6:51.001.57

Г. А. Щеглов

О ПРИМЕНЕНИИ ВОРТОННЫХ РАМОК В МЕТОДЕ ВИХРЕВЫХ ЧАСТИЦ

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

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

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

Для исследования пространственных течений вихревыми методами используются трехмерные вихревые элементы — вортоны: вихревая точка [3], вихревой сгусток (vortex blob) [4], вихревой отрезок [5]. Настоящая работа опирается на метод [6], в котором используется сферический вортон.

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

указанной расчетной схемы является одним из узких мест вихревых методов.

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

Постановка задачи. Пусть тело движется в невязкой жидкости с постоянной плотностью р = const, заполняющей все безграничное пространство. Рассматривается нестационарное трехмерное поле скоростей V (V, t). Движение жидкости описывается уравнением неразрывности div V = 0 и уравнением Эйлера, которое для решения задачи вихревым методом записывается в форме Ламба:

_ _ dt=(«V) * о

где Q = rotV — завихренность.

Граничные условия заданы на бесконечности в виде

lim V (V, t) = Vзо = const

r—^^o

и на поверхности тела S как

V (Vs,t) - Vs (rs,t)

ns = 0,

где ^ (г^, ¿) — скорость движения точки с радиус-вектором г3 на поверхности тела; — внешняя нормаль к поверхности тела в этой точке.

При дискретизации уравнений движения жидкости согласно методу вихревых частиц [6] область, занятая завихренностью (вихревой след), аппроксимируется множеством из N дискретных вихревых элементов — вортонов. Вортон в данном случае представляет собой сферу, заполненную завихренностью. Параметрами ¿-го вортона являются радиус-вектор центра сферы Гг, радиус сферы ег и вектор завихренности Шг. Скорость, индуцируемая вортоном в точке наблюдения с радиус-вектором Г), определяется зависимостью

иг (Г), ^ = К (бг, Го - Гг) X Шг (¿) ,

где

K (e,r)={ 4п |r| V е 0, Vr = 0,

з Л-) , Г =0;

., . 2 x5 + 5 x3

f(x) =-5.

2(x2 + 1)2

Если известны параметры всех вортонов вихревого следа, то скорость жидкости в точке наблюдения вычисляется по формуле

N

V (Т0,г) = Ко + ^ иг (Т0,г). (2)

1=1

Давление на заданном шаге расчета в точке пространства с радиус-вектором г0 может быть найдено с использованием аналога интеграла Коши-Лагранжа, рассмотренного в работе [7]:

п г) + V (Г0,1)2

р (го,Ь) = Рж + р—--р-

2 г 2

N

- р1Ьот (г) + р ^ V (тг,г) • иг (го, г), (3)

г=1

где 1Ьогп — вклад в давление от вортонов, рожденных на текущем шаге.

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

Движение вортонов описывается системой дифференциальных уравнений

ПТ - —»

= V(т1,г), г = 1,...,М. (4)

аг

Изменение векторов завихренности вортонов описывает система дифференциальных уравнений, получаемая из уравнения (1)

аг=<5(*, г), г=1,...,ж (5)

При движении вихревого следа для выполнения теорем Кельвина и Гельмгольца необходимо изменение радиуса вортонов. В работе [6] предлагается метод "расплющивания ядра". Он заключается в том, что на шаге расчета вортон удлиняется в направлении вектора завихренности. Этот процесс описывается дифференциальными уравнениями

dLi _ Lo i dt \шг\

düji

dt

dei ei (t) dL

(6) (7)

аг 2 ь0 аг'

где Ь0 — диаметр вортона в начале шага интегрирования, Ь — длина вортона в направлении вектора завихренности в конце шага интегрирования.

Удлиненный вортон в конце шага интегрирования заменяется сферой равного объема. Скорректированный вектор завихренности определяется по формуле

1 ((+ ш) М «С + ^ |( ^ .

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

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

Применение вортонных рамок. В настоящее время для описания процесса генерации завихренности в методе вихревых частиц используются вспомогательные дискретные особенности: поверхностные источники [6], вихревые плитки [8], замкнутые вихревые рамки, составленные из вихревых отрезков [9]. С их помощью на поверхности тела формируется расчетная схема для удовлетворения граничного условия непротекания. Таким образом, в ходе расчета возникает необходимость перехода от вспомогательных вортонов на теле к вортонам, моделирующим вихревой след. Указанный переход достаточно сложен, может приводить к потере точности и неустойчивости счета, а также требует дополнительных вычислительных ресурсов. Цель настоящей работы — построение расчетной схемы на теле, не требующей применения вспомогательных вортонов. Для этого развиваются идеи метода дискретных вихрей с замкнутыми вихревыми рамками [9].

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

В предлагаемой расчетной схеме вихревые отрезки, составляющие рамку, аппроксимируются сферическими вортонами. При этом образуется новый объект — вортонная рамка (рис. 1).

Параметрами ¿-го рожденного вортона являются: точка рождения с радиус-вектором Г®, радиус сферы е® и вектор завихренности Ш®. Расчетная схема строится таким образом, что Г®, е® и направление

Рис. 1. Вортонная рамка

вектора задаются геометрией обтекаемого тела, а неизвестной величиной, определяемой из условия непротекания, является только мо—* B

дуль вектора .

Пусть j-я m-угольная вихревая рамка задана радиус-векторами вершин fjs (s = 1,..., m), контрольной точки, в которой задано условие непротекания, k0j, циркуляцией Г и внешней нормалью к поверхности в контрольной точке nj. Тогда параметры вихревых элементов вортонной рамки определяются выражениями

B ' Ljs j £ B Ljs

rjsz = rjsz + —Trj's + 0 nj, £jsZ = ~ —T,

j = rjj; z = i,...,us.

где Ljs — длина, а fjs — единичный направляющий вектор s-го вихревого отрезка; us — число вортонов, аппроксимирующих отрезок.

В отличие от метода замкнутых вихревых рамок, где рассматриваются фиксированные линии отрыва вихревого следа, для вортон-ных рамок используется модель "потока завихренности" (vorticity flux) [10], в которой завихренность генерируется по всей поверхности обтекаемого тела и все вортоны каждой рамки после рождения пополняют вихревой след. Такая модель эффективна при решении задачи о колебаниях упругих тел вращения в жидкости, когда обтекаемая упругая деформируемая поверхность является гладкой и заранее невозможно указать области схода вихрей в поток. Для реализации "потока завихренности" в соотношения (8) вводится параметр 0 < S < min ef, который задает возвышение вортона над телом.

Как следует из соотношений (8) вортоны, рождающиеся над одной

f

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

U + J' us + 1' (8)

контрольных точках k0j

N

С

(9)

m us

где Wi (r0) = ^^ ^^ KK , r0 — rfsz) x Tls — вектор влияния вортон-

s=1 z=1

ной рамки; V ( k0j J — скорость жидкости в контрольной точке, вы-

числяемая по формуле (2); ^ — скорость тела в контрольной

точке, заданная уравнениями движения.

Для замкнутых тел матрица системы (9) получается вырожденной. В этом случае используется регуляризация, применяемая в методе дискретных вихрей [1, 9].

Алгоритм решения задачи. В основе алгоритма решения задачи лежит конвейер из пяти операций.

1. Построение расчетной схемы из п вортонных рамок на поверхности тела. На первом этапе поверхность тела разбивается на рамки. Определяются число сторон, координаты вершин, координаты контрольных точек и нормалей в них. На втором этапе вычисляются параметры вортонных рамок (8). На третьем этапе решается система алгебраических уравнений (9) и находятся неизвестные циркуляции

2. Пополнение вихревого следа новыми вортонами.

3. Восстановление поля давлений на поверхности тела по формуле (3) и вычисление гидродинамических нагрузок на данном шаге расчета.

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

5. Интегрирование уравнений эволюции вихревого следа. Все дифференциальные уравнения интегрируются численно методом первого порядка. Вортоны перемещаются в соответствии с (4), (5) и трансформируются согласно (6), (7). При перемещении вихревых частиц проверяется условие проницаемости. Вортоны, попавшие внутрь тела, возвращаются в поток, с использованием алгоритма, применяемого в методе дискретных вихрей [1].

Цикл вычислений продолжается до достижения конечного времени вычислений или до остановки вычислений пользователем программы.

Как показывают методические расчеты, для устойчивого счета достаточно применять вортонные рамки с параметром п3 = 1 (один

Г (j = 1,...,n).

вортон на стороне). Для снижения размерности задачи эффективным оказался прием объединения вортонов, принадлежащих соседним рамкам, после решения системы (9). В этом случае векторы завихренности вортонов складываются.

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

Результаты тестирования метода. Для тестирования метода была решена модельная задача с использованием безразмерных параметров. В вычислительном эксперименте изучалось обтекание жесткой неподвижной сферы радиусом Я = 0,5. В начальный момент времени среда (р = 1) приходит в движение со скоростью = { 1, 0, 0 }т.

Сфера разбита на п = 162 рамки (две десятиугольных рамки в полюсах и 160 четырехугольных рамок на остальной поверхности). Возвышение вортонов над телом принято равным 8 = 0,01.

Интегрирование системы дифференциальных уравнений эволюции завихренности проводилось с шагом = 0,04. Критерием останова расчета считалось установление коэффициента сопротивления. Установление считалось достигнутым в момент времени если на интервале — 50Д£; среднее квадратическое отклонение зависимости Сж(£) не превышало а = 0,05. В проведенном расчете критерий выполнялся уже на 150-м шаге расчета, однако для контроля устойчивости решения счет продолжался еще 20 шагов, в течение которых критерий не был нарушен. На рис. 2 показан вид вихревого следа, сформированного вортонами, число которых N = 13840, после 170 шагов расчета (£ = 6,8).

Оценка достоверности результатов проводилась путем сравнения с известными экспериментальными данными. На рис. 3 мгновенные линии тока, вычисленные на 170-м шаге, наложены на фотографию осредненного потока при обтекании шара для числа Рейнольдса Яе = 15000 [12]. Видно хорошее качественное совпадение картин течения.

Рис. 2. Вид вихревого следа, сформированного вортонами

Рис. 3. Вычисленные линии тока в сравнении с экспериментом [12]

• v Расчет ч / • / ° А ф

/ /

У* 90 V • \\ \ \ / 180 г • / \ Re= \ 270 =1,65-105 \ i \ / 360 • 7 ч

\ Re=4,35«105

Рис. 4. Распределение давления в сравнении с экспериментом [13]

На рис.4 показано распределение давления на поверхности для докритического и за-критического течения [13]. Были построены 60 эпюр p(^), вычисленных на каждом шаге по времени в интервале t Е [4,6; 6,8]. Для моделирования эффекта демпфирования пульсаций давления измерительной аппаратурой в численном эксперименте давление для каждой точки эпюры p(^) было взято средним на интервале t Е [4,6; 6,8]. Эти средние значения указаны на рис.4 точками.

На рис. 5 дан график изменения во времени осевой составляющей силы сопротивления в сравнении с результатом для Re = 15000 из работы [13]. Среднее значение, найденное по результатам вычислений, СХ = 0,57. Разница с экспериментальным значением СХэ = 0,49 составляет 16%.

Расчет, в котором было выполнено 170 полных циклов вычислительного конвейера, занял на персональном компьютере с процессором Pentium IV (2,6 ГГц) 16 часов. Для решения аналогичной нестационарной задачи методом конечных объемов потребовалось на том же компьютере 120 часов расчета.

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

I

V 1

1 1 1 1 1 2 3 4 t

Рис. 5. Вычисленные значения коэффициента сопротивления в сравнении с результатами эксперимента (штриховая линия) [13]

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

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

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

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

Работа поддержана грантом РФФИ № 06-08-00646.

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

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.Новиков Е. А. Обобщенная динамика трехмерных вихревых особенностей (вортонов) // Журнал эксп. и теор. физики. - 1983. - T. 84. - Вып. 3. - С. 975-981.

4. Winckelmans G. S., Leonard A. Contributions to vortex particle methods for the computation of Three-Dimensional Incompressible Unsteady Flows // Journal of comp. Physics, 109, 1993. - P. 247-273.

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

6. Ojima A., Kamemoto 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.

7. Д ы н н и к о в а Г. Я. Силы, действующие на тело при нестационарном вихревом отрывном обтекании идеальной несжимаемой жидкостью // Изв. РАН МЖГ. -2001. -№ 2. - С. 128-138.

8. Gharakhani A., Ghoniem A. F. Simulation of three dimensional internal flows by the random vortex and boundary element methods // ESAIM: Proceedings, Vol. 1, 1996.-P. 33-47.

9. Апаринов В. А., Дворак А. В. Метод дискретных вихрей с замкнутыми вихревыми рамками // Применение ЭВМ для исследования аэродинамических характеристик летательных аппаратов. - М.: Труды ВВИА им. проф. Н.Е. Жуковского, 1986. - Вып. 1313.

10. Chorin A. J. Vortex sheet approximation of boundary layers // Journal of computational physics. - 1978. - No. 27. - С. 428-442.

11. Альбом течений жидкости и газа: Пер. с англ. / Сост. М. Ван-Дайк. - М.: Мир, 1986. - 184 с.

12. Шлихтинг Г. Теория пограничного слоя. - М.: Наука, 1974. - 716 с.

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

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

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 21 publications in the field of dynamics and strength of spacecrafts.

ЖУРНАЛ "ВЕСТНИК МОСКОВСКОГО ГОСУДАРСТВЕННОГО ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА имени Н.Э. БАУМАНА"

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

Журнал "Вестник МГТУ имени Н.Э. Баумана" в соответствии с постановлением Высшей аттестационной комиссии Федерального агентства по образованию Российской Федерации включен в перечень периодических и научно-технических изданий, в которых рекомендуется публикация основных результатов диссертаций на соискание ученой степени доктора наук.

Подписку на журнал "Вестник МГТУ имени Н.Э. Баумана" можно оформить через агентство "Роспечать".

Подписывайтесь и публикуйтесь! Подписка по каталогу "Газеты, журналы" агентства "Роспечать"

Индекс Наименование серии Объем выпуска Подписная цена (руб.)

Полугодие 3 мес. 6 мес.

72781 "Машиностроение" 2 250 500

72783 "Приборостроение" 2 250 500

79982 "Естественные науки" 2 250 500

Адрес редакции журнала "Вестник МГТУ имени Н.Э. Баумана": 105005, Москва,

2-я Бауманская ул., д.5.

Тел.: (499) 263-62-60; (499) 263-67-98.

Факс: (495) 261-45-97.

E-mail: press@bmstu.ru

It*