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

Наука к Образование

МГТУ им. Н.Э. Баумана

Сетевое научное издание

Наука и Образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2014. № 12. С. 128-136.

Б01: 10.7463/0815.9328000

Представлена в редакцию: Исправлена:

© МГТУ им. Н.Э. Баумана

УДК 532.582.33

Вычисление распределения давления по поверхности тела вращения методом вихревых элементов

Дергачев С. А.

##.##.2014 ##.##.2014

1МГТУ им. Н.Э. Баумана, Москва, Россия

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

Ключевые слова: аэродинамика, верификация модели, вихревой метод, вихревой элемент, вортон, гидродинамика

Введение

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

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

На кафедре «Аэрокосмические системы» МГТУ им. Н.Э. Баумана развивается модификация метода вихревых элементов, где в качестве вихревого элемента используется симметричный вортон-отрезок [6]. Для удовлетворения граничных условий на поверхности используются замкнутые вортонные рамки [7]. Для расчета нестационарных нагрузок

на трехмерное тело при обтекании его дозвуковым потоком несжимаемой среды создана программа MVE3D которая. хорошо зарекомендовала себя при работе с телами, имеющими удлинение Ь / О < 6.

Целью данной отчета является проведение методических расчетов, направленных на выбор наилучших параметров метода вихревых элементов для моделирования распределения давления по поверхности затупленных тел вращения большого удлинения (Ь/О > 10). В таких телах, как показывают эксперименты [8], в зоне перехода между затуплением и цилиндрической частью возникает повышенное давление, которое требуется корректно моделировать для получения адекватных значений аэрогидродинамических характеристик..

1. Постановка задачи и метод решения

В неподвижной системе координат ОХУ2 рассматривается пространственное обтекание покоящегося абсолютно жесткого тела потоком несжимаемой среды, имеющим малую вязкость у << 1. Переход к безразмерным величинам производится так, что среда на бесконечном удалении от тела имеет плотность рх = 1,0, давление рш = 1,0, скорость

= (1,0;0,0;0,0}г. Обтекаемое тело представляет собой цилиндр с безразмерными параметрами: диаметр d = 0,2, длина Ь = 2,55. Передняя (наветренная) часть цилиндра имеет эллиптическое затупление с длиной большой полуоси Ь0 = 0,12. Тело расположено относительно системы координат ОХУ2 так, что угол атаки составляет а = 10°.

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

У-К = 0,

д¥

Г „Л

— + (у-УУ = УУ2К-V ^ V(тк) = 0, НшГ = К,Ншр = р,

(1)

где г - радиус-вектор точки в системе координат ОХУ2, г =| г |, гк - радиус-вектор

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

Допущение о несжимаемости среды позволяет определять поле скоростей по полю завихренности О. = Ух¥

V = К(О) + ^, (2)

где V(О) - поле скоростей, восстановленное с использованием закона Био-Савара.

г ^да

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

только как причину генерации завихренности ÙK в тонком пристеночном слое вблизи поверхности тела возникающей при удовлетворении граничного условия прилипания. Гипотеза Лайтхилла о потоке завихренности со всей поверхности обтекаемого тела [9] позволяет определить интенсивность генерируемой завихренности из условия непротекания поверхности тела.

V (r к ) • nK = (V (Q) + V (ÙK ) + Vx ) • nK = 0, (3)

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

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

DF=у

D (4)

DQ=(âv)V

Dt

с начальным условием Q = Q0, где r - радиус-вектор лагранжевой точки среды в

системе координат OXYZ.

Для решения (4) используется численный метод вихревых элементов. Поле завихренности представляется в виде суперпозиции NV элементарных полей - вихревых элементов (ВЭ)

Nv ( \

Q(r, t )«£По,. (г, 4 (t ), h (t ),r(t ))

i=1

Каждый ВЭ характеризуется тремя параметрами: маркером r0j (t), движущимся по траектории жидкой частицы, вектором h (t), и интенсивностью Г.. (t) . В качестве вихревого элемента используется симметричный вортон-отрезок для которого найдены аналитические выражения для скорости V(г, r0()Д- (t),r(t)) и ее градиента VV(r, r0()Д (t),T(t)), описанные подробно в [4]. Для исключения неограниченного роста скоростей и их производных при приближении к оси ВЭ, вводится радиус ВЭ s - трубки, внутри котороой индуцированные скорости и их производные убывают по линейному закону до нуля на оси ВЭ.

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

dr[

0г = V (4, t),

dt

^ = ^(тт)\ Иг, (г = 1,..., ^) (5)

dГ, л

—L = 0, dt

Начальными условиями для системы уравнений (5) являются параметры ВЭ в момент их рождения в процессе генерации завихренности. Генерация завихренности моделируется рождением N ВЭ из М вортонных рамок, образующих расчетную схему на поверхности тела, как описано в [5]. Вортонные рамки обеспечивают замкнутость и соленои-дальность вихревой системы вблизи поверхности. Каждая - -тая т. -угольная вихревая

рамка задана радиусами-векторами вершин (, = 1,...,т-), радиусом-вектором контрольной точки к^, внешней нормалью й^. к поверхности в контрольной точке и интенсивностью вихревого слоя у .

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

Мм=V}, (6)

где - матрица нормальных составляющих скосов потока, элементы которой равны

а у V] (к/ (здесь V,- (г) - вектор влияния ВЭ единичной интенсивности, лежащего

«=1

на ребре рамки, начинающемся в вершине ), {у} - вектор интенсивностей, {Уп} - вектор

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

{у} = [°У{УЯ}, (7)

После решения (6) происходит разделение рамок на ВЭ. При этом на ребрах соседних рамок рождается один, объединенный ВЭ. Параметры новых ВЭ, генерируемых в момент времени t = ^ выражаются через параметры вортонных рамок

4 = АГг (г, ), Иг = Лы (Г,- ),Ц. = АГ . (у, ), (. = 1,..., Np = 1,..., М) (8)

Размерность системы (5) при численном интегрировании увеличивается на каждом шаге из-за генерации новых ВЭ

N &+1) = ^ (tk) + Np

Для новых уравнений системы (5) начальными условиями служат параметры (8). Давление в контрольных точках рамок определяется при помощи аналога интеграла Коши-Лагранжа

ВЭ [10]. Из (9) видно, что чем дальше от тела находится ВЭ, тем меньший вклад он вносит в давление на поверхности тела. Поскольку наибольший вклад в давление вносят ВЭ, находящиеся вблизи тела, при решении задачи ВЭ, выходящие за пределы шара с центром в центре масс тела и радиусом Ь^аг, удаляются из расчетной схемы, что позволяет уменьшить рост размерности системы (5). Также рост размерности (5) ограничивается введением реструктуризации вихревого следа: соседние ВЭ объединяются.

Таким образом, численное решение задачи (1) на интервале времени 0 < / < /к сводится к выполнению вычислительного цикла с шагом по времени А/, содержащего N = ^ / А/ шагов. На каждом шаге выполняются следующие операций:

1. Вычисляется вектор {V}

2. По формуле (7) находится вектор {/}

3. Генерируется N ВЭ с параметрами (8) и увеличивается размерность (5)

4. Вычисляется распределение давления на теле по формуле (9)

5. Решается (5) методом первого порядка, определяются параметры ВЭ на следующем шаге

6. Проводится реструктуризация ВЭ для уменьшения размерности (5)

Для проведения расчетов использован алгоритм численного решения задачи, реализованный в программе МУЕЗБ.

На результаты, получаемые методом вихревых элементов, оказывают влияние основные параметры расчетной схемы: радиус ВЭ £ , определяющий гладкость поля скоростей, восстанавливаемого по закону Био-Савара в (2); шаг интегрирования (5), совпадающий с шагом расчета А/; форма, размеры и количество рамок, используемых для удовлетворения граничного условия на теле. Методические расчеты, выполненные ранее, показывают [], что первые два параметра зависят от характерного размера ребра рамки А1

Это показывает, что построение сетки из вортонных рамок на поверхности тела также важно, как построение сетки в расчетной области для сеточных методов. Поскольку от количества рамок М зависит скорость роста размерности системы (5), а, следовательно, и

где V (г) - скорость, индуцируемая / -тым ВЭ /в - интеграл, учитывающий генерацию

2. Результаты расчетов

— < £ < 2 А/

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

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

рамками в полюсах. Общее число рамок равно М = N • + 2.

Рис. 1. Набор рамок, моделирующих поверхность обтекаемого тела (N=85 и N0=20).

Получаемые в расчетах результаты сравнивались с результатами эксперимента, который был проведен в опытовом бассейне ЦНИИ имени академика А.Н.Крылова в 2004 году совместно с ОАО "ВПК "НПО Машиностроения". В ходе этого эксперимента были получены значения коэффициента давления Ср = 2(р-раУр^У^) в ряде точек на поверхности тела.

Также было проведено сравнение результатов, полученных методом вихревых элементов с результатами расчета распределения давления по поверхности тела, полученными сеточным методом конечных объемов (МКО) в программе SolidWorks FlowSimulation 2010. Моделирование было проведено в двух вариантах с разной дискретизацией сетки. Заметной разницы при увеличении дискретизации с 900000 ячеек до 1800000 ячеек не обнаружилось.

В ходе исследования рассматривались различные варианты дискретизации тела. Параметры лежали в пределах 45 < N < 119, 10 < NO < 30. Общее число рамок находилось в пределах 450<М <3570. Время расчета N = 3000 шагов составило 5<<35 часов.

Расчеты показали, что наименьшая дискретизация, позволившая получить результаты, близкие к данным эксперимента определялась параметрами N = 85, NO = 20. Это дает М = 1702 и А/ ~ 0,03. Отсюда радиус вихревого элемента принят равным е = А//2 = 0,015 , а шаг интегрирования Аt = 0,015. При меньшей дискретизации не удавалось вопроизве-сти пик давления в зоне перехода от затупления к цилиндрической части тела. При боль-

шеи дискретизации уменьшение шага интегрирования приводило к резкому росту размерности (5) и увеличению времени расчета поля скоростей по формуле (2) пропорционально

^3.

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

Ср от приведенной координаты Ь . Ь = —, где Ь - удаленность точки от сферического

затупления в продольном направлении, Ьт - длина тела.

Эпюры получены путем осреднения по заключительным шагам интегрирования Там же приведены эпюры, полученные методом конечного объема. Результаты полученные с помощью МУЕЗБ при достаточной дискретизации (N = 85, И0 = 20) обозначены

«МУЕЗБ», а при малой дискретизации ( N = 45, NO = 10 ) обозначены «МУЕЗБ-МБ».

а)

б)

Рис.2. Распределения давления на нижней (а) и верхней (б) образующих тела вращения

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

Заключение

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

Расчет удлиненных тел приводит к расчетным схемам с размерностью матрицы

системы алгебраических уравнений (6) M > 1700, что вызывает резкое увеличение времени счета и требует применения параллельных вычислений на большом количестве процессоров.

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

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

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

2. Калугин В.Т., Мордвинцев Г.Г., Попов В.М. Моделирование процессов обтекания и управления аэродинамическими характеристиками летательных аппаратов. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2011. - 527 с.

3. Alkemade. A.J.Q. On Vortex Atoms and Vortons: PhD Thesis. - Delft, (The Netherlands), 1994. - 209 p.

4. Cottet G.-H., Koumoutsakos P. Vortex Methods: Theory and Practice. - Cambridge: Cambridge University Press, 2000. - 320 p.

5. Katz J., Plotkin A. Low-speed aerodynamics. - Cambridge: Cambridge University Press, 2001. - 629 p.

6. Marchevsky I.K., Scheglov G.A. Symmetrical vortex fragmenton as a vortex element for incompressible 3D flow simulation (2011) Computational Fluid Dynamics 2010 - Proceedings of the 6th International Conference on Computational Fluid Dynamics, ICCFD 2010 PP. 897 - 898

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

8. Петров К.П. Аэродинамика транспортных космических систем. - М.: Эдиториал УРСС, 2000. - 366 с.

9. Lighthill M.J. Introduction. Boundary Layer Theory // Laminar Boundary Layers / Edited by J. Rosenhead. - New-York: Oxford University Press, 1963. - P. 54-61.

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

Science and Education of the Bauman MSTU, 2014, no. 12, pp. 128-136.

DOI: 10.7463/0815.9328000

Received: Revised:

##.##.2014 ##.##.2014

Science^Education

of the Bauman MSTU

I SS N 1994-0408 © Bauman Moscow State Technical Unversity

Calculation of Pressure Distribution at Rotary Body Surface with the Vortex Element Method

S.A. Dergachev

:Bauman Moscow State Technical University, Moscow, Russia

Keywords: aerodinamic, verification of model, vortex method, vortex element, vorton,

hydrodynamics

Vortex element method allows to simulate unsteady hydrodynamic processes in incompressible environment, taking into account the evolution of the vortex sheet, including taking into account the deformation or moving of the body or part of construction.

For the calculation of the hydrodynamic characteristics of the method based on vortex element software package was developed MVE3D. Vortex element (VE) in program is symmetrical Vorton-cut. For satisfying the boundary conditions at the surface used closed frame of vortons.

With this software system modeled incompressible flow around a cylindrical body protection elongation L / D = 13 with a front spherical blunt with the angle of attack of 10 We analyzed the distribution of the pressure coefficient on the body surface of the top and bottom forming.

The calculate results were compared with known Results of experiment.

Considered design schemes with different number of Vorton framework. Also varied radius of VE. Calculation make possible to establish the degree of sampling surface needed to produce close to experiment results. It has been shown that an adequate reproducing the pressure distribution in the transition region spherical cylindrical surface, on the windward side requires a high degree of sampling.

Based on these results Can be possible need to improve on the design scheme of body's surface, allowing more accurate to describe the flow vorticity in areas with abrupt changes of geometry streamlined body.

References

1. Trehmernoe otrivnoe obtekanie tel proizvolnoy formi [Three dimensions stalled flow around of body with different forms]/ Redactor S.M.Belocerkovskiy. - Moscow: TSAGI, 2000. -265p.

2. Kalugin V.T., Mordvincev G.G., Popov V.M. Modelirovanie processov obtekaniya i upravlenia aerodinamicheskimi harakteristikami letatelnih apparatov [Model of flow processes and control aerodynamics characteristics of fling bodies]-Moscow: Bauman MGTU Publ., 2011. - 527p.

3. Alkemade. A.J.Q. On Vortex Atoms and Vortons: PhD Thesis. - Delft, (The Netherlands), 1994. - 209 p.

4. Cottet G.-H., Koumoutsakos P. Vortex Methods: Theory and Practice. - Cambridge: Cambridge University Press, 2000. - 320 p.

5. Katz J., Plotkin A. Low-speed aerodynamics. - Cambridge: Cambridge University Press, 2001. - 629 p.

6. Marchevsky I.K., Scheglov G.A. Symmetrical vortex fragmenton as a vortex element for incompressible 3D flow simulation (2011) Computational Fluid Dynamics 2010 - Proceedings of the 6th International Conference on Computational Fluid Dynamics, ICCFD 2010 PP. 897 - 898

7. Scheglov G.A. O primenenii vortonnih ramok v metode vihrevih chastic [About using of vortex frame in vortex element method], Bauman MGTU Publ. Serial «Естественные науки». -2008. - №2. - С.104-113

8. Petrov K.P. Aerodinamika transportnih kosmicheskih sistem [Aerodinamic of transport space systems] Мoscow: Endotrial URSS, 2000. - 366p.

9. Lighthill M.J. Introduction. Boundary Layer Theory // Laminar Boundary Layers / Edited by J. Rosenhead. - New-York: Oxford University Press, 1963. - P. 54-61.

10. Andronov P.R., Guvernuk S.V., Dinnikova G. YA., Vihrevie metodi rascheta nestacionarnih gidrodinamicheskih nagruzok [Vortex methods of calculation non-stationary hydrodynamics forces]- МMoscow: Moscow unta, 2006. - 184 p.