CC BY
6
DOI: 10.15514/ISPRAS-2021-33(5)-15
Особенности реализации эффективного алгоритма параллельных вычислений для моделирования обледенения стреловидного крыла с профилем GLC-305
1'2 К.Б. Кошелев, ORCID: 0000-0002-7124-3945<koshelevkb@mail.ru> 2 А.В. Осипов, ORCID: 0000-0001-9223-4274 <a.osipov@ispras.ru> 2 С.В. Стрижак, ORCID: 0000-0001-5525-5180 <s.strijhak@ispras.ru>
1 Институт водных и экологических проблем СО РАН, 656038, Алтайский край, г. Барнаул, ул. Молодежная, д. 1 2 Институт системного программирования им. В.П. Иванникова РАН, 109004, Россия, г. Москва, ул. А. Солженицына, д. 25
Аннотация. В работе рассматривается возможность библиотеки ICELIB, разработанной в ИСП РАН, для моделирования процессов ледообразования на поверхности летательных аппаратов. В качестве тестового примера для сравнения точности моделирования физических процессов, возникающих при эксплуатации самолета, исследовалась поверхность стреловидного крыла с профилем GLC-305. В статье обсуждаются возможности эффективного алгоритма распараллеливания с использованием модели жидкой пленки, динамической сетки и геометрического метода биссектрис. Разработанная библиотека ICELIB это совокупность трех решателей. Первый решатель iceFoam1 предназначен для предварительной оценки зон обледенения поверхности фюзеляжа и крыла летательного аппарата. Изменением геометрической формы исследуемого тела пренебрегаем, толщина образования льда пренебрежимо мала. Данная версия решателя не имеет ограничений на количестве вычислительных ядер при распараллеливании. Вторая версия решателя iceDyMFoam2 предназначена для моделирования образования двух типов наледи гладкой ("Glaze ice"), так и рыхлой ("Rime ice") для которой зачастую форма льда принимает сложный и причудливый вид. Учитывается влияние изменения формы тела на процесс обледенения. Ограничения связаны с особенностями построением сетки вблизи пограничного слоя обтекаемого тела. Для движения передней и задней границ пленки используются разные алгоритмы, которые оптимизированы для своих случаев. Прирост производительности ограничен и достигается при фиксированном числе ядер. Третья версия решателя iceDyMFoam3 также позволяет учитывать влияние изменения поверхности твердого тела при образовании наледи на сам процесс обледенения. Для случая образования наледи гладкого типа последняя версия решателя пока уступает по своим возможностям второй при сложных формах поверхности льда. В третьей версии пока используется несколько упрощенный и более единообразный подход для расчета движения обеих границ пленки наледи. Проведена оценка результатов расчета с данными эксперимента M. Papadakis для различных профилей и стреловидного крыла для случая "Rime ice". Получено хорошее согласование с результатами эксперимента.
Ключевые слова: газокапельный; поток; крыло; расчет; решатель; сетка; обледенение; жидкая пленка; параллельные вычисления; координаты; узел; толщина льда; область; декомпозиция; осреднение; распределение по ядрам.
Для цитирования: Кошелев К.Б., Осипов А.В., Стрижак С.В. Особенности реализации эффективного алгоритма параллельных вычислений для моделирования обледенения стреловидного крыла с
профилем GLC-305. Труды ИСП РАН, том 33, вып. 5, 2021 г., стр. 249-258. DOI: 10.15514ASPRAS-2021-33(5)-15
Благодарости: Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 19-29-13016.
Features of the implementation of an efficient parallel computation algorithm for modeling the icing of a swept wing with a
GLC-305 airfoil
u K.B. Koshelev, ORCID: 0000-0002-7124-3945 <koshelevkb@mail.ru> 2A.V. Osipov, ORCID: 0000-0001-9223-4274 <a.osipov@ispras.ru> 2 S.V. Strijhak, ORCID: 0000-0001-5525-5180 <s.strijhak@ispras.ru> 1 Institute for Water and Environmental Problems SB RAS 1, Molodyoznaya St., Barnaul, 656038, Altai Krai, Russia 2 Ivannikov Institute for System Programming of the Russian Academy of Sciences, 25, Alexander Solzhenitsyn st., Moscow, 109004, Russia
Abstract. The paper considers the possibility of the ICELIB library, developed at ISP RAS, for modeling ice formation processes on the surface of aircraft. As a test example to compare the accuracy of modeling the physical processes arising during the operation of the aircraft, the surface of a swept wing with a GLC-305 profile was studied. The possibilities of an efficient parallelization algorithm using a liquid film model, a dynamic mesh, and the geometric method of bisectors are discussed. The developed ICELIB library is a collection of three solvers. The first solver iceFoam1 is intended for preliminary estimation of the icing zones of the fuselage surface and aircraft's swept wing. The change in the geometric shape of the investigated body is neglected, the thickness of ice formation is negligible. This version of the solver has no restrictions on the number of cores when parallelizing. The second version of solver iceDyMFoam2 is designed to simulate the formation of two types of ice, smooth ("Glaze ice") and loose ("Rime ice"), for which the shape of ice often takes on a complex and bizarre appearance. The effect of changing the shape of the body on the icing process is taken into account. The limitations are related to the peculiarities of the construction of the mesh near the boundary layer of the streamlined body. Different algorithms are used to move the front and back edges of the film, which are optimized for their cases. The performance gain is limited and is achieved with a fixed number of cores. The third version of solver iceDyMFoam3 also allows you to take into account the effect of changes in the surface of a solid during the formation of ice on the icing process itself. For the case of smooth ice formation, the latest version of the solver is still inferior in its capabilities to the second one with complex ice surface shapes. In the third version, a somewhat simplified and more uniform approach is still used to calculate the motion of both boundaries of the ice film. The estimation of the calculation results with the data of the experiment from M. Papadakis for various airfoils and swept wing for the case of "Rime ice" is carried out. Good agreement with the experimental results was obtained.
Keywords: gas-droplet flow; wing; simulation; solver; grid; ice accretion; liquid film; parallel simulation; coordinates; grid; thickness of ice; domain; decomposition; averaging; distribution on cores
For citation: Koshelev K.B., Osipov A.V. Strijhak S.V. Features of implementation an efficient parallel computation algorithm for modeling the ice accretion of a swept wing with a GLC-305 airfoil. Trudy ISP RAN/Proc. ISP RAS, vol. 33, issue 5, 2021, pp. 249-258 (in Russian). DOI: 10.15514/ISPRAS-2021-33(5)-15
Acknowledgements. The reported study was funded by RFBR, project number № 19-29-13016.
1. Введение
Вопросы изучения процессов обледенения являются актуальными в связи с использованием авиационной техники, различного оборудования, ветроэлектрических установок в северных широтах и в сложных климатических условиях на территории РФ. На поверхности элементов
самолета, вертолета, кабелей линий электропередач, лопастей ветроэлектрической установки (ВЭУ) могут образовываться различные формы льда, которые приводят к снижению аэродинамического качества и к снижению подъемной силы для обтекаемой поверхности. Для борьбы со льдом используют различные противобледенительные системы, но их не всегда устанавливают, например, на крыле самолета и лопасти ВЭУ. Поэтому важно иметь экспериментальные установки, технологии обработки данных летного эксперимента и расчетные эффективные коды для изучения процессов образования льда на поверхности при изменении физических параметров (скорость, температура, водность, диаметр капель, время в потоке капель) в широком диапазоне изменении исследуемых величин.
2. Математическая модель
В ИСП РАН был разработан тематический решатель iceFoam с использованием Эйлер-Лагранжевой модели для моделирования процесса обледенения [1]. В дальнейшем были разработаны три новые версии решателя iceFoam на основе Эйлерового подхода для газовой фазы и Лагранжевого подхода для моделирования капелек воды и модели пленки воды на поверхности профиля поверх слоя льда. Во всех трех версиях решателей используется модель пленки воды и слоя льда, которая связана с моделью по теории мелкой воды [1,2]. Первая версия решателя iceFoam1 предназначена для начальной оценки мест обледенения произвольного 3D профиля. Данный решатель полностью использует все возможности технологии параллельного программирования MPI для проведения эффективных вычислений средствами открытого пакета OpenFOAM v1912. Главным ограничением этого решателя является предположение, что толщина слоя льда достаточно мала и изменением поверхности профиля при обледенении можно пренебречь. Профиль со льдом и профиль без льда полагаются совпадающими. Для проведения декомпозиции пространственной области пленки используется разработанная утилита extrudeToFilmCellDist. Данная утилита позволяет обеспечить однообразие декомпозиции доменов в газовой области и в прилегающей пленке. Данная версия решателя не имеет ограничений на количестве ядер. С целью оценки масштабирования решателя выполнялись расчет на 12, 24, 48 вычислительных ядрах. С помощью данного решателя можно проводить быстрый расчет. Вторая версия решателя iceDyMFoam2 предназначена для моделирования обледенения как рыхлым льдом (Rime ice), так и гладким льдом (Glaze ice), для которого зачастую форма льда принимает сложный и причудливый вид. В этой версии решатель достаточно эффективно позволяет учитывать влияние изменения поверхности профиля при обледенении на сам процесс обледенения. Данный решатель требует использования структурированной сетки вблизи поверхности профиля и наличия только одного пространственного домена для пленки и примыкающей к ней эйлеровой сетки. Опыт вычислений показал, что это требование ограничивает максимальное количество ядер, не более 12. При дальнейшем же увеличении количества используемых ядер эффективность параллельных растет слабо или даже падает при 24 ядрах. Для упрощения декомпозиции пространственной области газа разработана утилита addTwoLayersToO. Данная утилита позволяет задать расчетный домен на нулевом процессоре для пленки и для двух прилегающих слоев эйлеровой сетки газовой фазы потока. Третья версия решателя iceDyMFoam3 также позволяет учитывать влияние изменения поверхности профиля при обледенении на сам процесс обледенения. Как и первая версия, этот решатель полностью использует все возможности технологии параллельной библиотеки MPI для проведения эффективных параллельных вычислений. Параллельная реализация границы движения льда потребовала использовании низкоуровневых функций на языке программирования C++ открытого пакета OpenFOAM. В этом случае будут обязательно узлы, которые попадут на границу нескольких доменов.
В случае использования подхода с динамической сеткой потребовалось собирать информацию об узлах, которые хранятся в нескольких доменах или на разных расчетных вычислительных ядрах. Для таких узлов после локальных вычислений внутри одного домена выполнялась операция по определению нового положения узла. Далее делалось осреднение по всем доменам, участвующих в расчетах данного узла. Осредненные значения положение узлов, которые записаны были больше чем в одном домене, обратно распределялись по всем необходимых вычислительным ядрам. Таким образом, обеспечивалось совпадения движущихся узлов на разных доменах.
Для случая образования наледи гладкого типа третья версия решателя iceDyMFoam3 пока уступает по своим возможностям второй версии решателя при сложных формах поверхности льда, так как во второй версии для движения передней и задней границ пленки используются разные алгоритмы, которые оптимизированы для своих случаев. В третьей же версии пока используется несколько упрощенный и более единообразный подход для расчета движения обеих границ. Для проведения декомпозиции пространственной области пленки используется разработанная утилита extrudeToFilmCellDist.
Все три версии решателя iceFoam и дополнительные утилиты для работы с расчетной сеткой составляют основу тематической библиотеки ICELIB.
Для верификации решателей нами была проведена серия расчетов на вычислительных кластерах ИСП РАН, НИЦ "Курчатовский Институт" для различных профилей (NACA 0012, GLC-305, General Aviation, Business Jet, Commercial Transport) и 3D стреловидного крыла. Расчеты проводились для различных значений диаметров капель, температуры, скорости набегающего потока, водности, времени обледенения в набегающем потоке. Было проведено сравнение расчета с результатами эксперимента по толщине и форме льда [3,4]. Ошибка не превышала более 10%.
3. Особенности реализации параллельного алгоритма
Согласно концепции пакета OpenFOAM, узлы на границах расчетных доменов дублируются в каждом из этих доменов. В случае изменения координат этих узлов необходимо обеспечить их совпадение на разных вычислительных ядрах сервера. Имеется значение допустимого рассогласования координат таких узлов. Обычно это значение 0.0001 м. Для расчета изменения границы льда геометрическим методом биссектрис необходима информация о положении всех граней, которым принадлежит узел. Однако, по умолчанию, такая информация доступна только в рамках текущего домена. Фактически на границе доменов вместо процедуры интерполяции используется процедура экстраполяция. При этом, как показано на рис. 1, перемещение совпадающих узлов из разных доменов происходит неоднообразно.
correctly
Рис. 1. Распределение по доменам и направление вектора роста льда Fig. 1. Domain distribution and direction of the ice growth vector
Поэтому для устранения этого недостатка была разработана и запрограммирована
следующая процедура:
252
1. Предварительный шаг. Выявляются номера всех узлов, которые принадлежат более чем одному домену. Этот шаг выполняется один раз перед расчетом, так как топология сетки не меняется.
2. На каждом шаге по времени при расчете движения границы льда:
2.1. Стандартно выполняется процедура расчета движения льда методом биссектрис.
2.2. Для каждого узла, принадлежащего более чем одному домену, выполняется осреднение координат по всем ядрам процессоров с использованием функции "reduce".
2.3. Осредненное значение координаты присваивается соответствующим узлам также для каждого ядра.
4. Постановка задачи
Для расчетного примера рассматривалась стреловидное крыло с профилем GLC-305 из работы [5]. На рисунке 2 показана геометрия стреловидного крыла с длиной L=1.52 метра и шириной корневой хорды c=0.64 метра, на рисунке 3 показан профиль крыла с координатами точек в безразмерном виде. Расчетная область представляла собой параллелепипед с размерами: длина 7 метров (от -2 до 5 метров, крыло начинается в точке с координатами 0, 0, 0), высота 3 метра (от -1.5 до 1.5 метров), ширина 3 метра (от 0 до 3 метра). Для крыла также заданы углы стреловидности и крутки.
На всех внешних границах, кроме входного сечения, задавалось граничное условие "zero gradient", но поверхности крыла - условие "no slip". В одном парселе получилось 350000 капелек, а всего парселей было больше миллиона. Сферические частицы вводились за 0.5 метра от передней корневой кромки крыла в сечении по вертикали от -0.19 до -0.01 метра, чтобы накрыть крыло полностью под углом в 6 градусов, по ширине 1.6 метра. Ранее подобные численные расчеты проводились с использованием различных собственных кодов для аналогичного стреловидного крыла в работах [6-8].
Рис. 2. 3D модель стреловидного крыла Fig. 2. 3D swept wing model
S 1 L lr
Рис.3 Модель профиля GLC-305 Fig.3.Profile model GLC-305
Расчетная неструктурированная адаптивная сетка была подготовлена с помощью утилиты snappyHexMesh в составе ОрепРОАМ v19.12 и содержала 1 355 049 узлов. Сетка включала в себя две области: первая для пристеночной подобласти пленки с высотой в одну расчетную ячейку, вторая для внешней области для моделирования газокапельного потока. Характерный размер ячейки около крыла был порядка 4 мм. На рис. 4 изображена сетка с 4 уровнями адаптации.
Рис. 4. Сетка с 4 уровнями адаптации Fig. 4. Mesh with 4 levels of adaptation
5. Результаты расчета
Расчеты проводились для случая "Rime ice" для значения скорости V=90 м/c, угла атаки в 6 градусов, температура воздуха Tair = 261.87 K, LWC=0.51 г/м3, диаметр частиц d=14.5 микрон, время обледенения tau=5 минут. Расчетный шаг задавался равным 6-10"5 секунд. Один вариант считался на вычислительном кластере HPC4 НИЦ "Курчатовский Институт" 33 часа на 24 вычислительных ядрах.
a) b)
Рис. 5. Поле скорости и траектории частиц вокруг стреловидного крыла Fig. 5. Velocity field and particle trajectories around a swept wing
a) b)
Рис. 6. Форма льда на крыле в разных сечениях A Fig. 6. The shape of the ice on the wing in different sections A
a) b)
Рис. 7. Форма льда на крыле в разных сечениях B Fig. 7. The shape of the ice on the wing in different sections B
a) b)
Рис. 8. Форма льда на крыле в разных сечениях C Fig. 8. The shape of the ice on the wing in different sections C
Рис. 9. Форма льда на нижней поверхности крыле Fig. 9. Ice shape on the underside of the wing
На рис. 5 показана вертикальная скорость капель воды с видом сверху и снизу. Скорость относительно хорды крыла, которая сама наклонена под углом атаки 6 градусов. На рис. 6а, 7а, 8а передняя кромка крыла до и после обледенения примерно в тех же местах, что и в [5]. На рис. 6b, 7b, 8b весь профиль до и после обледенения в тех же местах. Видно, что к краю крыла обледенение снизу увеличивается.
На рис. 9 представлена форма поверхности льда снизу на крыле. Видно 3D положение поверхности льда. Этот рисунок демонстрирует способность решателя моделирования существенно трехмерных процессов обледенения.
В работе [5] приведены осредненные только 2D графики поверхности льда. При этом можно отметить, что результаты данной работы близки по максимальной толщине льда (примерно 3 мм) к данным эксперимента, по форме льда результаты также качественно совпадают.
6. Заключение
В данной работе представлено описание трех версий решателя iceFoam c поддержкой параллельных вычислений. Решена задача с моделированием обледенения для стреловидного крыла с профилем GLC-305. Для расчетов использовалась версия решателя iceDyMFoam3. Проведено сравнение с результатами эксперимента работы [5]. Получено хорошее совпадение по толщине и форме льда.
Необходима дальнейшая работа по улучшению распараллеливания алгоритма решателя iceDyMFoam3 для расчетов на многоузловом вычислительном кластере c использованием технологии MPI [9,10].
Список литературы / References
[1]. Кошелев К.Б., Мельникова В.Г., Стрижак С.В. Разработка решателя iceFoam для моделирования процесса обледенения. Труды ИСП РАН, том 32, вып. 4, 2020 г., стр. 217-234 / Koshelev K.B., Melnikova V.G. Strijhak S.V. Development of iceFom solver for modeling ice accretion. Trudy ISP RAN/Proc. ISP RAS, vol. 32, issue 4, 2020. pp. 217-234 (in Russian). DOI: 10.15514/ISPRAS-2020-32(4)-16.
[2]. Bourgault Y., Beaugendre H., Habashi W. Development of a shallow-water icing model in FENSAP-ICE. Journal of Aircraft, vol. 37, no. 4, 2000, pp. 640-646.
[3]. Wright W., Rutkowski A. Validation results for LEWICE 2.0. Technical report CR-1999-208690, NASA, 1999, 674 p.
[4]. Addy, Harold E. Ice Accretions and Icing Effects for Modern Airfoils. Technical report TP-2000-210031, NASA, 2000
[5]. Papadakis M., Yeong H-W. et al. Experimental investigation of ice accretion effects on a swept wing. Technical report PB2005-110681, NASA, 2005, 205 p.
[6]. Zocca M., Gori G., and Guardone A. Blockage and Three-Dimensional Effects in Wind-Tunnel Testing of Ice Accretion over Wings. Journal of Aircraftm, vol. 54, no. 1, 2016, pp. 759-767.
[7]. Pena D., Hoarau Y., Laurendeau E. A single step ice accretion model using Level-Set method. Journal of Fluids and Structures, vol. 65, 2016, pp. 278-294.
[8]. Chang S., Tang H. et al. Three-Dimensional Modelling and Simulation of the Ice Accretion Process on Aircraft Wings. International Journal of Astronautics and Aeronautical Engineering, vol. 3, issue 2, 2018, pp. 1-25.
[9]. Гергель В.П. Высокопроизводительные вычисления для многопроцессорных многоядерных систем. М., Изд-во МГУ, 2010 г., 544 стр. / Gergel V.P. High performance computing for multiprocessor multicore systems. M., Publishing house of Moscow State University, 2010, 544 p. (in Russian).
[10]. Гергель В.П., Сысоев А.В. Высокопроизводительные параллельные вычисления. 100 заданий для расширенного лабораторного практикума. М., ФИЗМАТЛИТ, 2018 г., 248 стр. / Gergel V.P., Sysoev A.V. High performance parallel computing. 100 tasks for an extended laboratory practice. M., FIZMATLIT, 2018, 248 p.
Информация об авторах / Information about authors
Константин Борисович КОШЕЛЕВ - кандидат физико-математических наук, доцент, старший научный сотрудник. Сфера научных интересов: вычислительная гидродинамика, гидрология, геоинформатика.
Konstantin Borisovich KOSHELEV is a candidate of physical and mathematical sciences, associate professor, senior researcher. Research interests: computational fluid dynamics, hydrology, geoinformatics.
Андрей Владимирович ОСИПОВ - инженер. Сфера научных интересов: вычислительная гидродинамика, метод контрольного объема, подвижные сетки, лагранжев подход.
Andrei Vladimirovich OSIPOV - an engineer. Research interests: computational fluid dynamics,
finite volume method, dynamic meshes, particles.
256
Сергей Владимирович СТРИЖАК - кандидат технических наук, ведущий инженер. Сфера научных интересов: вычислительная гидродинамика, многофазные течения, турбулентность, ветроэнергетика, параллельные вычисления.
Sergei Vladimirovich STRIJHAK - a candidate of technical sciences, leading engineer. Research interests: computational fluid dynamics, multiphase flows, turbulence, wind power, parallel computing.
