CC BY
111
Исследование применения солверов пакета ОрепГОЛМ для моделирования сверхзвукового обтекания конуса
Бондарев А.Е., Кувшинников А.Е., ИПМ им. М.В. Келдыша РАН
bond@keldysh.ru
Аннотация
Приводятся численные результаты сравнения точности для ряда солверов пакета ОрепРОЛЫ. Проведено сравнение для моделирования сверхзвукового обтекания конуса невязким сжимаемым газом при нулевом угле атаки. Результаты, полученные с помощью различных солверов ОрепРОЛЫ, сравниваются с известным численным решением задачи при вариации угла конуса и скорости потока. Данная работа является частью проекта, направленного на создание надежной вычислительной технологии для моделирования потоков вокруг удлиненных тел вращения (УТВ).
1 Введение
Методы численного моделирования задач аэродинамики в настоящее время находят широкое применение на практике. Существует большое число численных моделей и основанных на этих моделях программных комплексов для решения задач такого вида. При этом становится актуальным вопрос о применимости подобных моделей, как правило, состоящих из большого числа параметров, подбор которых объясняется чисто эмпирическим путем или вовсе никак не объясняется. Эта проблема актуальна в различных областях научных исследований.
Исследованию этой проблемы в некотором частном случае, а именно для задач обтекания удлиненных тел вращения, и посвящена данная статья.
Задача обтекания УТВ привлекательна тем, что для некоторых постановок существуют как аналитическое решение, так и достаточное количество экспериментальных данных.
Требуется вычислительная технология, которая была бы своего рода эталоном для решения задач обтекания УТВ и помогала регулировать настраиваемые параметры как численных методов, так и моделей турбулентности в различных программных пакетах. В этом качестве было решено воссо-
здать на уровне современных высокопроизводительных вычислительных средств вычислительную технологию, разработанную ранее в ИПМ им. М.В. Келдыша А.Е. Бондаревым и В. А. Черкашиным под руководством А.В. Забродина. Данная вычислительная технология позволяла в конце 80-х - начале 90-х годов надежно проводить массовые промышленные расчеты УТВ практического назначения. Погрешность коэффициентов аэродинамического сопротивления не превышала 2-3 процентов по сравнению с экспериментальными результатами. Суть данной технологии заключалась в том, что коэффициент аэродинамического сопротивления УТВ Сх, необходимый для расчетов баллистики исследуемого объекта, рассматривался как сумма трех составляющих: Cp - сопротивления УТВ при невязком обтекании, Cf - коэффициента сопротивления трения и Сд - донного сопротивления. Подобный подход был широко распространен в задачах массового промышленного анализа аэродинамических свойств УТВ и показывал себя весьма эффективным [1].
Часть этой вычислительной технологии уже реализована. Так, например, для определения коэффициента трения реализована вычислительная методика [2-3], основанная на приближенной полуэмпирической модели, объединяющей результаты экспериментальных исследований Л.В. Козлова [4-5] и метод эффективной длины, предложенный В.С. Авдуевским [6-7].
Для расчета аэродинамических характеристик УТВ при невязком обтекании было предложено использовать программный пакет OpenFOAM (Open Source Field Operation And Manipulation CFD Toolbox) [8]. Это свободно распространяемый программный продукт, написанный на языке C++. OpenFOAM активно используется в промышленности и в науке. OpenFOAM содержит в себе ряд солверов, обладающих различными вычислительными свойствами.
Стоит пояснить, что в OpenFOAM существуют как заранее созданные разработчи-
ками солверы [9-11], так и возможность создавать собственные солверы, что использовали авторы из ИСП РАН [12-13].
В качестве точного решения традиционно использовались табличные решения [14]. Решения, представленные в [14], получены с помощью конечно-разностных методов для обтекания гладких тел потоком невязкого газа в широком диапазоне чисел Маха и углов полураствора конуса с вариацией угла атаки. Решения представлены в виде таблиц, обладают высокой точностью и на протяжении многих лет используются в качестве точного решения при анализе вычислительных свойств того или иного численного метода.
Следует отметить, что подобные сравнения солверов проводились в работах [1517]. Однако эти сравнения проводились на других примерах и не дают четких рекомендаций по выбору солвера для рассматриваемого класса задач.
2 Постановка задачи
Постановка задачи представлена в полном соответствии с работой [14], где рассматриваются результаты невязкого обтекания конусов с различными углами полураствора на различных числах Маха.
Исследуется обтекание УТВ равномерным сверхзвуковым потоком идеального газа под нулевым углом атаки а = 0° при числе Маха М ¥ = 2. Исследуемое тело является конусом с углом полураствора Ь = 10°- 25° с шагом 5° . Условия набегающего потока на входе обозначаются индексом «1», а на выходе — индексом X, так как решение является автомодельным и зависит от безразмерной переменной. Схема течения представлена на рисунке 1.
Для расчета используется система уравнений Эйлера. Система дополняется уравнением состояния идеального газа.
Рис. 1. Схема течения
3 Выбор солверов OpenFOAM
Для сравнения из программного пакета OpenFOAM были выбраны 4 солвера.
A) rhoCentralFoam — основан на цен-трально-противопотоковой схеме, которая является комбинацией центрально-разностной и противопотоковой схем [910]. Суть центрально-противопотоковых схем состоит в специальном выборе контрольного объема, содержащего области двух типов: вокруг граничных точек - первый тип; вокруг центральной точки - второй тип. Границы контрольных объемов первого типа определяются при помощи локальных скоростей распространения возмущений. Преимущество указанных схем состоит в том, что, применяя соответствующую технику уменьшения численной вязкости, можно добиться хорошей разрешимости и для разрывных решений - ударных волн в газовой динамике - и для решений, где основную роль играют вязкие явления.
Б) sonicFoam — основан на алгоритме PISO (Pressure Implicit with Splitting of Operator) [11]. Основная идея метода PISO заключается в том, что для расчета давления используются два разностных уравнения для поправки поля давления, полученных из дискретных аналогов уравнений моментов и неразрывности. Такой подход связан с тем, что скорректированные первой поправкой скорости могут не удовлетворять уравнению неразрывности, поэтому вводится второй корректор, который позволяет вычислить скорости и давления, удовлетворяющие линеаризованным уравнениям количества движения и неразрывности.
B) rhoPimpleFoam — основан на алгоритме PIMPLE, который является комбинацией алгоритмов PISO и SIMPLE (Semi-Implicit Method for Pressure-Linked
Equations) [18]. К алгоритму PISO добавляется внешний цикл, благодаря которому метод становится итерационным и позволяет считать с числом Куранта, большим 1.
Г) pisoCentralFoam — комбинация цен-трально-противопотоковой схемы с алгоритмом PISO [12].
Расчеты для всех солверов проводились с помощью программного пакета OpenFOAM версии 2.3.0.
4
4.1
Организация расчетов и полученные результаты
Построение сетки, начальные и граничные условия
На рисунке 2 представлена расчетная область. Постановка граничных условий представлена в таблице 1. На верхней границе, обозначенной в таблице «top», задается условие нулевого градиента для газодинамических функций, обозначаемое в таблице «zeroGradient». Такие же условия задаются на правой границе, обозначаемой «outlet». На левой границе, обозначаемой «inlet», заданы параметры набегающего потока: давление P = 101325 Па, температура T = 300 K, скорость U = 694.5 м/с. На границе конуса «cone» для давления и температуры задается условие нулевого градиента, для скорости задается условие «slip», соответствующее условию непротекания для уравнений Эйлера. Для моделирования осесимметричной геометрии в пакете OpenFoam для передней «front» и задней «back» границ используется специальное условие «wedge». Для оси «axis» в пакете OpenFoam также вводится специальное граничное условие «empty». Это условие задается в случаях, когда вычисления в заданном направлении не проводятся.
Чтобы оценить влияние размера ячейки сетки на точность, расчеты для U=2M проводились на трех сетках, обозначенных как coarse, fine, finest. Число ячеек сетки coarse 3000, fine - 12000, finest - 48000. Так как размер области остается неизменным, то при увеличении числа ячеек сетки, размер ячеек уменьшается.
Рис. 2. Расчетная область
Табл. 1. Граничные условия
P T U
inlet 101325 300 Mach 2 -Mach 5
outle t zeroGradien t zeroGradient zeroGradient
axis empty empty empty
top zeroGradien t zeroGradient zeroGradient
cone zeroGradien t zeroGradient slip
front wedge wedge wedge
back wedge wedge wedge
4.2 Параметры солверов
В пакете OpenFOAM существует два варианта аппроксимации дифференциальных операторов: непосредственно в коде решателя или с использованием файлов конфигурации fVSchemes и fvSolution. Для того, чтобы сравнение было правильным, мы использовали одинаковые параметры там, где это было возможно.
В файле fVSchemes: ddtSchemes - Euler, gradSchemes - Gauss linear, divSchemes -Gauss linear, laplacianSchemes - Gauss linear corrected, interpolationSchemes- vanLeer.
В файле fVSolution: solver - smoothSolver, smoother symGaussSeidel, tolerance - 1e-09, nCorrectors - 2, nNonOrthogonalCorrectors -1.
4.3 Расчет осесимметричного течения
На рисунках 3 и 4 представлены поля установившегося течения для давления и плотности при использовании солвера rhoCentralFoam. Эти рисунки свидетельствуют о том, что в результате установления получена качественная картина тече-
ния, соответствующая решениям [14].
известным
Рис. 3. Поле давления для установившегося течения
Табл. 2. Отклонение от точного решения для сетки coarse
rCF
0.009062
0.043725
0.024054
0.018327
pCF
0.008929
0.050789
0.027705
0.021848 0.028965
sF
0.008366
0.050932
0.033429
rPF
0.010155
0.060268
0.037406
0.033199
Табл. 3. Отклонение от точного решения для сетки fine
Uy
rCF
0.006268
0.029656
0.016989
0.012834
pCF
0.006482
0.034403
0.019515
0.015182 0.019085
sF
0.005809
0.033814
0.022465
rPF
0.007588
0.043562
0.026656
0.022994
Табл. 4. Отклонение от точного решения для сетки finest
ux Uy
rCF
0.004372
0.019862
0.011611
0.008715
pCF
0.004441
0.022855
0.013269
0.010282
sF
0.004057
0.023113
0.015143
0.012684
rPF
0.005526
0.030994
0.018803
0.015810
U
U
p
p
U
p
p
Рис. 4. Поле плотности для установившегося течения
В таблицах (2-8) приведены результаты расчетов в форме аналога нормы Ь2:
exact У m.
exact У m.
где ут —компоненты скорости их, Цу, давление р и плотность р в ячейке, Ут— объем ячейки для угла полураствора в = 1035° с шагом 5° и чисел Маха М=2-5. Жирным выделены минимальные значения. Символ "х" в таблицах означает, что при данной скорости и данном угле полураствора солвер становился неустойчивым. Здесь значения утехаС получены интерполированием табличных значений из [14] на ячейки сетки. Следует отметить, что авторы таблиц [14] указывают на допустимость интерполяции по всем параметрам и значениям таблиц.
Далее для солверов будем использовать сокращенные обозначения: гСР
(ЛоСепйаШоат), рСР (pisoCentralFoam), sF ^отсРоат), гРР (ЛоРтр1еРоат).
В таблицах (2-4) приведены результаты отклонения от точного решения для всех исследуемых солверов для случая в = 20° и числа Маха M=2. Таблицы содержат данные для всех газодинамических функций. Приведены результаты для трех сеток coarse, fine и finest соответственно. Результаты иллюстрируют сходимость численного решения.
Таблицы (5-8) представляют собой результаты для давления в форме аналога нормы L2, представленные при вариации угла полураствора от 10° до 35° и вариации числа Маха набегающего потока от 2 до 5.
Табл. 5. Отклонения от точного решения, U=2M
Угол полураствора
10 15
20 25 30 35
rCF
0.00609 0
0.01265 4
0.01662 3
0.01867 8
0.02069 5
0.03248 6
pCF
0.00697 3
0.01444 6
0.01935 3
0.02094 8
0.02313 0
0.03865 8
sF
0.01015
0.01964 6
0.02228
0.02077 9
0.02561 4
0.07484 9
rPF
0.01034 1
0.02064 5
0.02495 1
0.02542 6
0.02326 7
0.04317 9
2
m
Табл. 6. Отклонение от точного решения, и=3Ы
Угол полураствора гСТ pCF sF
10 0.01530 9 0.01953 7 0.02715 2 0.02717 7
15 0.02460 8 0.03004 1 0.04781 3 0.04144 4
20 0.03044 0 0.03585 8 0.07056 4 0.04576 0
25 0.03248 6 0.03865 8 0.07484 9 0.04317 9
30 0.03404 0 0.04060 3 0.07740 8 0.04000 6
35 0.02633 4 0.02982 1 0.04485 3 0.02707 7
0,025-| 0,02 0,015 0,01 0,005
0-М
¿=1
гСР рСР sF
гРР
Табл. 7. Отклонение от точного решения, и=4Ы
Угол полураствора гСТ pCF sF
10 0.02825 4 0.03525 1 0.05813 3 0.04933 4
15 0.04022 9 0.04649 4 0.10617 2 0.06538 4
20 0.04615 9 0.05268 7 0.12670 1 0.07064 9
25 0.04584 9 0.05191 2 0.13493 2 0.06278 5
30 0.04077 5 0.05061 9 0.10912 5 X
35 0.03427 7 0.04229 6 0.06966 8 X
Рис. 5. Отклонение от точного решения для давления при M=2, в = 20
На рисунке 6 представлено изменение отклонения от точного решения в аналоге нормы Ь2 для давления для всех солверов в зависимости от угла полураствора конуса при фиксированном числе Маха Ж=2. Наименьшее отклонение от точного решения показывает солвер гЬоСеПхаШоаш, наибольшее отклонение при увеличении угла полураствора показывает солвер 8ошсРоаш.
Табл. 8. Отклонение от точного решения, и=5Ы
Угол полураствора гСТ pCF sF
10 0.05083 4 0.05513 3 0.10671 0 0.07582 9
15 0.06006 9 0.06329 3 0.15988 0 0.09048 9
20 0.06017 4 0.06467 5 0.17566 6 X
25 0.05990 0 0.06328 4 0.17520 5 X
30 0.05597 5 0.06263 7 0.13020 1 X
35 0.04328 8 0.05273 7 0.09000 6 X
Рисунок 5 представляет диаграмму отклонения от точного решения в аналоге нормы L2 для давления для всех использованных солверов на примере задачи обтекания конуса с углом полураствора в = 20° при числе Маха M=2. Наименьшее отклонение от точного решения показывает сол-вер гЬоСеПгаШоаш, наибольшее отклонение показывает солвер гЬоРтр1еРоаш.
угол в
Рис. 6. Отклонение от точного решения для давления в зависимостиот угла конуса для всех солверов, M = 2
На рисунке 7 представлена зависимость отклонения от точного решения в аналоге нормы Ь2 для давления для солвера гЬоСеПхаШоаш при вариации угла полураствора и начальной скорости. Увеличение числа Маха набегающего потока оказывает наибольшее влияние на увеличение отклонения численного результата от точного решения.
0,08 0,06 0,04 0,02 0
5 M
10 15 20 25 30 35 ß
Рис. 7. Отклонение от точного решения
для давления в зависимости от угла полураствора конуса и скорости для солвера ЛоСейгаШоат
5 Заключение
Для задачи обтекания конуса сверхзвуковым потоком идеального газа под нулевым углом атаки было проведено сравнение четырех солверов программного пакета ОрепБоат с точным решением, полученным из таблиц [14]. Результаты расчетов были представлены в целях сравнительного анализа в форме аналога нормы Ь2.
Согласно полученным результатам сол-вер АоСеп^аШоат обладает минимальной нормой погрешности поля почти во всех случаях. Единственным недостатком ЛоСейгаШоат является появление осцил-ляций у поверхности в головной части конуса. Солвер pisoCentrlFoam находится на втором месте по точности, однако при использовании этого солвера появления ос-цилляций не наблюдается.
Для числа Маха М=2 солвер soпicFoam более точно вычисляет компоненты скорости. Но для числа Маха М=3 и выше в результатах, полученных при помощи этого солвера наблюдаются осцилляции на фронте ударной волны. Это негативно сказывается на норме погрешности. Солвер rhoPimpleFoam не работает для числа Маха большего М=4 при углах полураствора больше в = 20°. Кроме того, по точности он почти никогда не превосходит солверы rhoCentralFoam и pisoCentrlFoam. Таким образом, можно утверждать, что солверы rhoCentralFoam и pisoCentrlFoam обеспечивают лучшую точность для интересующего класса задач и могут быть использованы при построении вычислительной технологии расчетов обтекания удлиненных тел вращения.
Проведенное методическое исследование может служить основой для выбора солвера программного пакета OpenFoam при расчете невязкого сверхзвукового обтекания УТВ, а также может быть полезным для разработчиков программного контента OpenFoam.
В дальнейшем предполагается провести аналогичное исследование для более высоких чисел Маха M >5. Также предполагается провести аналогичное сравнение солве-ров для задачи обтекания конуса под углом атаки при вариации угла атаки, угла полураствора и числа Маха набегающего потока. Следует также заметить, что в сравнительный анализ могут быть дополнительно включены результаты для новых солверов, полученных от разработчиков.
Благодарности
Данная работа выполнена при поддержке Российского фонда фундаментальных исследований (проекты 16-01-00553а, 17-01-00444а, 18-31-00320мол-а).
Список литературы
[1] Красильщиков А.П., Гурьяшкин Л.П. Экспериментальные исследования тел вращения в гиперзвуковых потоках. М: ФИЗМАТЛИТ, 2007. 208 с.
[2] Приближенный подход к оценке сопротивления трения на телах вращения в вязком потоке / С.В.Андреев [и др.] // Препринты ИПМ им. М.В.Келдыша. 2014. № 102. 12 с. URL:
http://library.keldysh.ru/preprint. asp?id=2014-102.
[3] Bondarev A.E., Nesterenko E.A. Approximate method for estimation of friction forces for ax-isymmetric bodies in viscous flows // Mathemat-ica Montisnigri. Vol. XXXI. 2014. P.54-63.
[4] Козлов Л.В. Экспериментальное исследование поверхностного трения на плоской пластине в сверхзвуковом потоке при наличии теплообмена // Изв. АН СССР. Механика и машиностроение. 1963. №2. С.11-20.
[5] Козлов Л.В. Экспериментальное определение закона теплообмена для турбулентного пограничного слоя в сверхзвуковом потоке // Исследование теплообмена в потоках жидкости и газа. М.: Машиностроение. 1965. С. 91-109.
[6] Основы теплопередачи в авиационной и ракетной технике [Текст] / В.С. Авдуевский, Ю.И. Данилов, В.К. Кошкин и др.; Под общ. ред. проф. В. К. Кошкина. Москва: Оборон-гиз, 1960. 389 с.
[7] Авдуевский В.С. Метод расчета пространственного турбулентного пограничного слоя в сжимаемом газе // Изв. АН СССР. Механика и машиностроение. 1962. №4. С.3-13.
[8] OpenFOAM, URL: http://www.openfoam.org. (дата обращения 01.02.2018).
[9] Kurganov A., Tadmor E. New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations. // J. Com-put. Phys. 2000. Vol. 160. P. 241-282, doi:10.1006/jcph.2000.6459.
[10] Implementation of semi-discrete, non-staggered central schemes in a colocated, polyhedral, finite volume framework, for high-speed viscous flows / Christopher J. Greenshields, Henry G. Wellerr, Luca Gasparini, Jason M. Reese // Int. J. Numer. Meth. Fluids. 2010, Vol. 63. Issue 1. P. 1-21, doi:10.1002/fld.2069.
[11] Issa R. Solution of the implicit discretized fluid flow equations by operator splitting. // J. Comput. Phys. 1986. Vol. 62. Issue 1. P. 40-65, doi:10.1016/0021 -9991 (86)90099-9.
[12] Kraposhin M., Bovtrikova A., Strijhak S. Adaptation of Kurganov-Tadmor numerical scheme for applying in combination with the PISO method in numerical simulation of flows in a wide range of Mach numbers. // Procedia Computer Science. 2015. Vol. 66, P. 43-52. doi:10.1016/j.procs.2015.11.007.
[13] Github: United collection of hybrid Central solvers
https://github.com/unicfdlab/hybridCentralSolve rs (дата обращения 01.02.2018).
[14] Пространственное обтекание гладких тел идеальным газом / К. П. Бабенко, Г. П. Воскресенский, А. Н. Любимов, В. В. Русанов. М.: Наука, 1964.
[15] Численное моделирование сверхзвукового обтекания клина с применением свободного открытого программного кода OpenFOAM. / А. Я. Карвацкий, П. В. Пулинец, Т. В. Лазарев, А. Ю. Педченко // Космiчна наука i тех-нолопя. 2015. Т. 21, № 2. С. 47-52.
[16] Gutierrez L. F. M., Tamagno J. P., Elaskar S. A. High speed flow simulation using Open-FOAM. // Mecanica Computacional. 2012. Vol XXXI. P. 2939-2959.
[17] Lorenzon D., Elaskar S. A. Simulacion de flujos supersonicos bidimensionales y axialmen-te simetricos con OpenFOAM. // Revista de la
Facultad de Ciencias Exactas, Fisicas y Naturales, 2015. Vol. 2. no. 2. P. 65-76. URL: http://revistas.unc .edu. ar/index. php/F CEF yN/article/view/11061
[18] Patankar S. V., Spalding D. B. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. // International Journal of Heat and Mass Transfer. 1972. Vol. 15, Issue 10. P. 1787-1806. doi: 10.1016/0017-9310(72)90054-3
