CC BY
36
Математическое моделирование
эволюции завихренности при пространственном обтекании тел методом вихревых петель
С.А. Дергачев <sadergachev@mail.ru> МГТУ им. Н.Э. Баумана, 105005, г. Москва, ул. 2-я Бауманская, д.5. стр.1
Аннотация. Моделирование гидродинамических явлений, связанных с обтеканием подвижных деформируемых тел является актуальной инженерной задачей для различных технических приложений. В данной статье описан метод вихревых петель, позволяющий моделировать обтекание тел без необходимости предварительного задания линий схода завихренности. Приведено верифицирование методики на экспериментальных данных обтекания крыла конечного размаха и сферы. Показана точность, достаточная для применения программы и методики в инженерных приложениях. Сделаны выводы о возможности перехода к подвижным и деформируемым телам. Программа реализована на языке С++ с использованием библиотеки MPI для организации пересылки данных при параллельных вычислениях.
Ключевые слова: гидродинамика; нестационарные нагрузки; вихревые петли; вихри; метод вихревых петель
DOI: 10.15514/ISPRAS-2018-30(1 )-14
Для цитирования: Дергачев С.А. Математическое моделирование эволюции завихренности при пространственном обтекании тел методом вихревых петель. Труды ИСП РАН, том 30, вып. 1, 2018 г., стр. 215-226. DOI: 10.15514/ISPRAS-2018-30(1)-14
1. Введение
Вычисление нестационарных гидродинамических нагрузок на подвижные деформируемые трехмерные объекты в безграничном объеме несжимаемой среды является актуальной задачей в различных инженерных приложениях. Существуют сеточные и бессеточные подходы к вычислению этого вида нагрузок. При условии значительных перемещений и деформаций тела нестационарные нагрузки с точки зрения быстродействия расчетов удобнее вычислять бессеточными методами вычислительной гидродинамики. В случае отрывного обтекания тела несжимаемой средой нестационарные нагрузки во многом определяются процессами вихреобразования. В вихревых методах
вычислительной гидродинамики непосредственно моделируется возникновение и эволюция завихренности, а поле скоростей и давление вычисляются по известному полю завихренности. Вихревые бессеточные методы различаются по способам моделирования завихренности. Наиболее известен метод дискретных вихрей, получивший развитие в научной школе С.М. Белоцерковского [1]. При моделировании пространственных течений завихренность заключена в замкнутых вихревых рамках, состоящих из наборов отрезков. Рамки в следе соединены в тонкие вихревые пелены. Узлы вихревых рамок перемещаются по полю скорости. Поверхность тела также моделируется замкнутыми вихревыми рамками. Рамки, моделирующие поверхность тела перемещаются в соответствии с перемещениями или деформациями тела. Рамки формируются при сходе завихренности с заданных линий на поверхности тела. Линии схода назначаются заранее перед проведением расчета, что затрудняет моделирование гладких деформируемых тел. Кроме того, как показывает практика расчетов, невозможность разрыва и перезамыкания вихревых пелен может служить причиной вычислительной неустойчивости в случае взаимодействия нескольких пелен или пелены с обтекаемой поверхностью.
В методах вихревых элементов завихренность моделируется отдельными, не связанными между собой, частицами - вортонами, компактными носителями завихренности. Основная часть завихренности в них сконцентрирована в ограниченной области: в точке, в эллипсоиде или на отрезке. При расчете обтекания тел обычно используется гипотеза Лайтхилла-Чорина, согласно которой вортоны генерируются по всей поверхности обтекаемого тела. Такой подход позволяет эффективно моделировать отрывы потока с гладких деформируемых тел. Однако, в силу теорем Гельмгольца во всем пространстве вокруг вортона распределена еще некоторая часть завихренности, называемая дополнительной. Дополнительная завихренность обеспечивает соленоидальность поля вихревого элемента. Ее интенсивность стремится к нулю на бесконечном удалении от компактной области. Однако результаты расчетов показывают, что дополнительная завихренность может приводить к неправильному вычислению нагрузок на цилиндрические тела [2]. Известно, что в природе завихренность имеет свойство образовывать замкнутые вихревые структуры: нити, кольца, петли. Данное свойство формализовано в теоремах Гельмгольца. Моделированию эволюции вихревых нитей, колец и петель посвящено значительное количество исследований [3,4]. Замкнутая непересекающаяся ломаная линия - вихревая петля, составленная из вортонов-отрезков, может быть использована в качестве вихревого «суперэлемента». Если вортоны, моделирующие такие элементы, имеют одинаковую интенсивность и связанны между собой дополнительная завихренность отсутствует, поскольку вся завихренность заключена в вихревых петлях.
Применение вихревых петель для расчета нестационарных гидродинамических нагрузок в инженерных задачах потребовало разработки
эмпирических моделей их эволюции, генерации на поверхности обтекаемого тела и вычисления на основе петель поля давления и нагрузок на конструкцию [5]. Целью данной работы является верификация разработанных моделей, алгоритмов и программы по известным экспериментальным данным, а также определение допустимых значений параметров расчета.
2. Постановка задачи
Моделируется течение несжимаемой среды, для которой эффекты вязкости учитываются только как причина рождения завихренности на поверхности тела и как причина перезамыкания вихревых петель.
Используются два уравнения гидродинамики: уравнение неразрывности и уравнение сохранения импульса (уравнение Гельмгольца):
(Цу 7 = 0; ^ = го£(7 х П);
где 7 (г, £) - нестационарное трехмерное поле скоростей; г- радиус-вектор точки в неподвижной системе координат; П = го£7 - завихренность. Заданы граничные условия отсутствия возмущений на бесконечности и условие прилипания на обтекаемой поверхности. Начальные условия соответствуют режиму бесциркуляционного обтекания тела. Требуется в течение заданного промежутка времени определить распределение давления по поверхности тела и вычислить действующие на него нагрузки.
В рамках вихревого метода уравнение неразрывности удовлетворяется тождественно, поскольку поле скорости восстанавливается по полю завихренности в соответствии с законом Био-Савара. Уравнение сохранения импульса заменяется системой обыкновенных дифференциальных уравнений, описывающей в лагранжевой постановке эволюцию завихренности, которая моделируется системой из Ыр вихревых петель, каждая из которых состоит из Ыр вортонов-отрезков. Все вихревые петли в расчёте имеют одинаковую интенсивность Г. Данная величина задается в составе исходных данных. Вклад в поле скоростей от одного вортона-отрезка единичной интенсивности определяется формулой:
Р /В _ 1 В (11 - А^ В2 • А^ ^
^ (Г)=^ рр[1В ¡ВТ
В = В X АВ , В = В-, Р2 = В- Р„+1, / = 1,...,к =
Для устранения сингулярности в расчетах вводится радиус линейного сглаживания скорости - радиус вихревого элемента е. Величины Ыр и Ыр изменяются в процессе расчета из-за рождения новых петель на обтекаемом теле, а также из-за изменения длины и перезамыкания петель. Блок-схема алгоритма решения задачи приведена на рис. 1.
Рис. 1. Схема организации работы программы. Жирным выделены распараллеленные
процедуры
Fig. 1. Scheme of the organization of the program. Parallel procedures are marked in bold.
Поверхность обтекаемого тела моделируется набором из N плоских панелей. Каждая панель представляет собой рамку из вортонов-отрезков, а в центре панели расположена контрольная точка. На каждом расчетном шаге из
условия обеспечения непротекания в контрольных точках (равенства нулю нормальной составляющей скорости) из решения СЛАУ находятся интенсивности рамок.
Далее распределение интенсивности панелей на поверхности тела сглаживается. По полученному скалярному полю, начиная от полюса с максимальной на данном шаге интенсивностью, осуществляется построение линий уровня с шагом Г. Вдоль этих линий уровня формируются новые петли. После создания петли сбрасываются в поток - отодвигаются от тела на некоторую величину Д.
Сформированные таким образом петли подвергаются нормализации: элементы петель попавшие внутрь тела из-за погрешностей интегрирования заменяются на элементы петель проложенные по поверхности на удалении Д, для обеспечения устойчивости расчета осуществляется устранение сильных изломов (т.н. шпилек) и проводится перераспределение узлов вдоль петель для того чтобы обеспечить длину отрезков равную заданному параметру a. Алгоритмы нормализации петель сформированы по эмпирическим законам, путём верификации расчетов с экспериментальными данными [5]. Величина р определяет минимальный угол между отрезками. Если угол между отрезками меньше, то данное место считается шпилькой и узел переносится ближе для уменьшения угла до допустимой величины. В расчетах получено оптимальное значение 160°.
Давление в контрольных точках определяются аналогом интеграла Коши-Лагранжа [6], адаптированным под метод вихревых петель:
_ ®-А-Ь/» - -У £/ЦЛ-1
р(Г, () =
V
2 2 4пА1 £ "" 1 ,
где А1 - телесный угол под которым видна панель из точки вычисления давления.
На следующем этапе расчёта определяются скорости в узлах вихревых петель. Численное интегрирование движения узлов по траекториям жидких частиц проводится методом Эйлера первого порядка в заданным шагом по времени dt. Для компенсации погрешностей интегрирования после перемещения узлов снова проводится коррекция фрагментов петель, пересекающих тело и нормализация петель.
В конце шага определяются близко расположенные участки петель. При близком расположении (ближе, чем величина ц) противоположно направленных участков петель осуществляется их перезамыкание. В результате петли могут объединяться или разделяться на части. Процедура перезамыкания сопровождается нормализацией петель и производится в цикле до тех пор, пока не будут обработаны все петли.
Описанный алгоритм реализован в виде программы MVortexLoops на языке С++ с использованием библиотеки MPI для распараллеливания вычислений. На разработанную программу получено свидетельство о государственной регистрации №2017616752. Возможность параллельных вычислений была реализована для наиболее трудоёмких процедур: определения скоростей; нормализации петель; коррекции элементов петель, пересекающих тело; перезамыкания петель.
Решение СЛАУ не распараллеливалось, поскольку перемещение и деформация тел на данном этапе не моделировалось и обращения матрицы системы проводилось лишь на первом шаге расчёта обтекания тел.
3.Результаты моделирования
На первом этапе верификации производилось моделирование обтекания крыла конечного размаха (длина крыла 0,5, хорды =0,1), имеющего профиль NACA-2217 и удлинение 5,0 потоком с единичной скоростью V(r, t) = {1,0,0} при различных углах атаки.
Поверхность модели разбита на четырехугольные и треугольные панели. Общее число панелей поверхности тела равно N = 5212. Разбиение торца профиля на панели приведено на рис. 2. Видно, что задняя кромка затуплена, для уменьшения взаимного влияния панелей на верхней и нижней поверхностях профиля, а дискретизация профиля в передней части является достаточно грубой, что связано с необходимостью соблюдать одинаковый размер ячеек для всей модели. Верхняя и нижняя поверхности разбиты на прямоугольные рамки. При этом ширина рамок соответствует стороне граничных треугольников на торцах крыла.
Рис. 2. Разбиение боковой поверхности профиля на панели Fig. 2. Wing lateral surface partition
В расчетах радиус вихревого элемента взят равным е=0,002, шаг по времени At=0,00219, длина отрезка a=0,001, интенсивность петель Г=0.00526, параметры алгоритма расчета эволюции петель Д = 0.0001, р = 160°, /л, = 0.04.
На рис. 3 приведены фазы самоорганизации за крылом вихревых жгутов из вихревых петель (красные линии): рождение и унос разгонного вихря и формирование за крылом П-образной системы Прандтля. Также на рисунке видна присоединенная к крылу завихренность в виде системы вихревых петель, а также отдельные вихревые петли, формирующие тонкую вихревую пелену за острой кромкой. Таким образом в расчетах с использованием 220
вихревых петель получена качественно верная картина формирования вихревого следа.
Рис. 3. Развитие вихревого следа за крылом Fig. 3. Vortex wake evolution behind the wing
На рис. 4 приведено сравнение рассчитанных в скоростной системе координат аэродинамических коэффициентов с экспериментальными данными [7]. Коэффициенты сил и момента вычисляются качественно верно, ошибка в определении подъемной силы и момента не превышают 8%, что допустимо для проведения инженерных расчётов. Относительная ошибка лобового сопротивления велика в силу малости величины лобового сопротивления профиля, однако абсолютное значение погрешности также допустимо для инженерных расчетов. Причина завышенного значения лобового сопротивления и заниженного значения подъемной силы, возможно, заключена в грубости модели профиля, примененного для расчетов. Здесь требуются дополнительные методические исследования. На втором этапе верификации было проведено моделирование обтекания сферы единичного диаметра. Поверхность сферы разбита на N=1917 четырехугольных панели. Вектор скорости набегающего потока был единичным 7(г, £) = {1,0,0}. Всего было выполнено около 520 расчетов. При этом варьировались параметры расчетной схемы: шаг по времени dt, радиус сглаживания скорости вблизи
петли е, дистанция перезамыкания петель ц, длина отрезка а, дистанция от поверхности тела до узлов петель при создании Д, интенсивность петель Г.
Мг,
-0т4—
Угол ä аки и. °
О -: 5 1 ) 5 2
■ —ьз-
Рис. 4. Аэродинамические коэффициенты. Вычисленные значения - точки, данные
экспериментов - линии [7] Fig. 4. Aerodynamics coefficient. Calculated value - points, experimental data - lines [7]
Результаты расчетов сравнивались с результатами распределения коэффициента давления, экспериментально полученными для 4 значений числа Рейнольдса вблизи кризиса сопротивления сферы, как показано на рис.5.
При расчетах обтекания сферы использовался кластер 10P МСЦ РАН. Каждый запуск программы проводился на 16 ядрах, а длительность счета соответствовала 960 минутам. Периоды вычисленного переходного режима (количество сделанных шагов) во всех расчетах различны и зависят от конкретных расчетных параметров.
Критерием качества итеративного подбора параметров расчетной схемы являлось среднеквадратичное отклонение точек №1,2,3,4 (см. Рис.5) в плоскости OXY значений коэффициента давления, полученных при расчете, от экспериментально полученных значений при соответствующих числах Рейнольдса. В результате для каждого режима обтекания получены наилучшие сочетания параметров, представленные в табл. 1. В таблице первый столбец соответствует номеру расчета, далее шесть параметров расчета, T - длительность переходного режима, y1...y4 - среднеквадратичные
отклонения результатов расчета от данных эксперимента. Распределения коэффициента давления для найденных сочетаний параметров показаны на рис. 6.
о* 20" 4П- SDn 8DJ ГОО" 120" МО' 100° 180° Ч>
Рис. 5. Распределение коэффициента давления вдоль образующей шара 1-Re=157200, cx=0,471; 2-Re=251300, cx=0,313; 3-Re=298500, cx=0,151; 4-Re=424500, cx=0,143 [8] Fig. 5. The distribution of the pressure coefficient along the generatrix of the sphere 1-Re=157200, cx=0,471; 2- Re=251300, cx=0,313; 3- Re=298500, cx=0,151; 4-Re=424500, cx=0,143 [8].
Результаты (см табл. 1) показывают, что разработанная методика позволяет моделировать гидродинамические нагрузки при больших числах Рейнольдса, когда влияние сил вязкости мало по сравнению с влиянием сил инерции в жидкости. При этом шаг по времени примерно соответствует отношению среднего линейного размера панели к скорости невозмущенного набегающего потока.
Данные рекомендации по выбору расчетных параметров соответствуют рекомендациям при расчёте методом вихревых элементов. При меньших значениях числа Рейнольдса количественно схожие результаты получаются при увеличении шага по времени, увеличении радиуса сглаживания и увеличении дистанции отодвигания петель от поверхности тела при рождении.
Табл. 1 Параметры и результаты моделирования обтекания сферы Table 1 Parameters and simulation results sphere in the flow
NN £ dt a Д Г T yi y2 y3 y4
482 0.0087 0.23 0.102 0.06 0.042 0.01 395 0.049 0.148 0.293 0.343
521 0.0005 0.10 0.052 0.12 0.016 0.01 195 0.1539 0.0631 0.168 0.214
550 0.0035 0.027 0.02 0.08 0.00125 0.02 21.6 0.287 0.150 0.050 0.076
309 0.0035 0.025 0.02 0.08 0.00125 0.02 28.2 0.3477 0.2109 0.090 0.046
NN=482 Re= 157200
NN=521 Re=251300
Рис. 6. Сравнение распределения коэффициента давления, полученного в расчете(точки) и экспериментально (линии) Fig. 6. Comparison of the pressure coefficient distribution obtained in the calculation (points) and experimentally (lines)
Следует учесть, что в разработанной методики влияние вязкости на взаимодействие между петлями не учитывается. Для корректного моделирования влияния вязкости на процессы рождения и взаимодействия вихревых петель методика должна быть доработана, например, с использованием модели обмена интенсивностями между вортонами. Однако, поскольку непосредственной целью разработанной методики является моделирование гидродинамических нагрузок при больших числах Рейнольдса (около 7 106), методика может быть использована без доработок.
4. Заключение
Проведенная верификация разработанных моделей, алгоритмов и программы MVortexLoops по известным экспериментальным данным показала, что использование замкнутых вихревых петель в качестве вихревого «суперэлемента» позволяет качественно верно моделировать нестационарное
отрывное обтекание крыльев и гладких тел. При этом количественные значения аэродинамических коэффициентов и распределения давления могут быть найдены с приемлемой для инженерных приложений погрешностью при соответствующем выборе параметров эмпирических моделей. Применение разработанной модификации метода вихревых элементов позволяет, с одной стороны, не задавать до начала расчета линии отрыва вихревых пелен, и, с другой стороны, позволяет исключить вредное влияние дополнительной завихренности вортонов-отрезков. Дальнейшее развитие будет направлено на моделирование гидродинамических нагрузок на подвижные деформируемые тела.
Работа выполнена с использованием средств Межведомственного суперкомпьютерного центра Российской академии наук.
Список литературы
[1]. Трехмерное отрывное обтекание тел произвольной формы. Под ред. С.М. Белоцерковского. М.: ЦАГИ, 2000, 265 с.
[2]. Щеглов Г.А. Модификация метода вихревых элементов для расчета гидродинамических характеристик гладких тел. Вестник МГТУ им Н.Э. Баумана. Сер. Машиностроение, 2009, №2, с. 26-35
[3]. Chorin, A.J. Hairpin removal in vortex interactions II, Jour. of Comput. Phys., 1993. №107, . p. 1-9,
[4]. WeibmannS., Pinkall U., Filament-based smoke with vortex shedding and variational reconnection. ACM Transactions on Graphics, 2010, Vol. 29, No. 4, Article 115. DOI: 10.1145/1778765.1778852
[5]. Дергачев С.А., Щеглов Г.А.. Моделирование обтекания тел методом вихревых элементов с использованием замкнутых вихревых петель. Научный вестник МГТУ ГА, 2016, № 223(1), с.19-25 DOI: 10.26467/2079-0619-2016--223-19-27
[6]. Андронов П.Р., Гувернюк С.В., Дынникова Г.Я. Вихревые методы расчета нестационарных гидродинамических нагрузок. М.: Изд-во Моск. ун-та, 2006, 184 с.
[7]. Б.А. Ушаков. Атлас аэродинамических профилей крыльев. Издание БНТ НКАП при ЦАГИ, 1940, 84 с.
[8]. Авиация: Энциклопедия. Гл. ред. Г.П. Свищёв М.: Большая Российская энциклопедия, 1994. 736 с.
Mathematical simulation of vorticity evolution in the case of spatial flow around bodies by the method of vortex loops
S.A. Dergachev <sadergachev@mail.ru> Bauman Moscow State Technical University, 5/1, 2-nd Baumanskaya st., Moscow. Russia, 105005
Abstract. Simulation of hydrodynamic phenomena associated with flow around movable deformable bodies is an actual engineering task of various technical problems. This article describes the method of vortex loops, which allows to simulate the flow around bodies
without the need to preset the lines of vorticity retreat. Verification of the technique on the experimental data on the flow past the wing and the sphere is given. The accuracy is shown sufficient for application of the program and methodology in engineering applications. Conclusions are made about the possibility of transition to mobile and deformable bodies. The program is implemented in C ++ with the use of the MPI library to organize the transfer of data in parallel calculations.
Keywords: hydrodynamics; non-stationary loads; vortex loops; vortices; method of vortex loops.
DOI: 10.15514/ISPRAS-2018-30(1)-14
For citation: Dergachev S.A. Mathematical simulation of vorticity evolution in the case of spatial flow around bodies by the method of vortex loops. Trudy ISP RAN/Proc. ISP RAS, vol. 30, issue 1, 2018, pp. 215-226 (in Russian). DOI: 10.15514/ISPRAS-2018-30(1)-14
References
[1]. Three-dimensional detached flow around bodies of arbitrary shape. S.M. Belotserkovsky Ed.. Moscow: TsAGI, 2000. - 265 p. (in Russian)
[2]. Shcheglov G.A. Modification of the vortex element method for calculating the hydrodynamic characteristics of smooth bodies. Vestnik MGTU im N.E. Baumana. Ser. Mashinostroenie [Herald of the Bauman Moscow State Technical University. Ser. Mechanical engineering]. 2009. № 2. p. 26-35 (in Russian)
[3]. Chorin, A.J. Hairpin removal in vortex interactions II. Jour. of Comput. Phys., 1993, №107, P. 1-9.
[4]. WeibmannS., Pinkall U., Filament-based smoke with vortex shedding and variational reconnection. ACM Transactions on Graphics, 2010. Vol. 29, No. 4, Article 115
[5]. Dergachev S.A., Scheglov G.A. Simulation of flow around the bodies by means of the method of vortex elements with use of closed vortex loops. Nauchny vestnik MGTU GA [Scientific Herald of the MGTU GA], 2016. No. 223(1), p.19-25. (in Russian). DOI: 10.26467/2079-0619-2016--223-19-27
[6]. Andronov P.R., Gouvernyuk S.V., Dynnikova G.Ya. Vortex methods for calculating unsteady hydrodynamic loads. Moscow: Izd-vo Mosk. University, 2006, 184 p. (in Russian)
[7]. Ushakov B.A. Atlas of airfoils. Edition of the BNT NCAP at TsAGI, 1940, 84 p. (in Russian)
[8]. Aviation: Encyclopedia. G.P. Svishchev Ed. Moscow: The Great Russian Encyclopedia, 1994, 736 p. (in Russian)
