МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

2019 год, том 31, номер 10, стр. 117-129



ВЫЧИСЛИТЕЛЬНЫЙ ЭКСПЕРИМЕНТ В ЗАДАЧЕ

СВЕРХЗВУКОВОГО ОБТЕКАНИЯ ЗАТУПЛЕННОГО ТЕЛА

С ХВОСТОВЫМ РАСШИРЕНИЕМ

2019 г.

И.А. Широков¹, Т.Г. Елизарова²

¹Московский государственный университет им. М.В. Ломоносова, факуль-

тет ВМК

ivanshirokov@inbox.ru

²Институт прикладной математики им. М.В. Келдыша РАН telizar@mail.ru

Работа поддержана грантом РФФИ 19-01-00262.

DOI: 10.1134/S0234087919100095

Представлены результаты прямого численного моделирования сверхзвукового об-

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

полученных данных с результатами лабораторного эксперимента для чисел Маха 3

и 4 и углов атаки в 10 и 20 градусов. Результаты получены с применением квазига-

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

Ключевые слова: сверхзвуковое обтекание, квазигазодинамические уравнения, не-

структурированные сетки.

COMPUTATIONAL EXPERIMENT IN THE PROBLEM

OF SUPERSONIC FLOW AROUND A BLUNT BODY WITH TAIL EXPANSION

I.A. Shirokov¹, T.G. Elizarova²

¹Moscow State University, Faculty of Computational Mathematics and Cybernetics

²Keldysh Institute of Applied Mathematics, RAS, Moscow

The results of direct numerical simulation of supersonic flow around a cylindrical blunt

body with tail expansion are presented. The data obtained are compared with the results

of laboratory experiments for the Mach numbers 3 and 4 and the angles of attack of 10

and 20 degrees. The results are obtained by using the quasi-gas-dynamic equations.

Keywords: supersonic flows, quasi-gas-dynamic equations, unstructured spatial grids.

1. Введение

В работе приводятся результаты прямого численного моделирования

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

трехмерной постановке с формированием вблизи ее поверхности отрывных

вихревых зон, которые могут вызывать генерацию акустических колебаний.

118

И.А. Широков, Т.Г. Елизарова

Для выбранной модели имеются данные лабораторных измерений и резуль-

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

сматриваться как тестовый вариант для оценки точности численных подхо-

дов. Расчет течения производится с использованием квазигазодинамических

(КГД) уравнений без применения традиционных для подобных задач вы-

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

было показано, что численный алгоритм, основанный на КГД уравнениях,

позволяет моделировать турбулентные течения при небольших числах Рей-

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

В представленных далее расчетах модели турбулентности также не приме-

няются.

Целью работы является изучение возможности применения КГД алго-

ритма, реализованного на неструктрурированной тетраэдральной расчетной

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

ний вязкого газа, в том числе при использовании не очень хороших сеток.

Модель представляет собой цилиндрическое затупленное тело с рас-

ширением в хвостовой части. Общий вид и размеры модели приведены на

рис.1 [1].

Рис.1. Общий вид и размеры модели.

2. Постановка задачи и газодинамические параметры

В стендовом эксперименте, описанном в [1], модель расположена в на-

бегающем потоке воздуха, при этом угол атаки принимает значения 10^ и

20^ . Рассматриваются числа Маха Ma 3 и Ma  4. Число Рейнольдса, от-

несенное к длине модели 0.127 м, равняется

Re  610

(при Ma  3) и Re 

 7.610

(при Ma  4). Отметим, что соответствующие числа Рейнольдса,

отнесенные к 1м, равны

4.710

(при Ma 3) и

610

(при

Ma  4). Параметры набегающего потока имеют следующие значения: газо-

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

119

вая постоянная R  287Дж/(кг·К), показатель адиабаты   1.4, число Пран-

дтля Pr  14 /19 , показатель межмолекулярного взаимодействия   0.74 .

3. Математическая модель

Система квазигазодинамических (КГД) уравнений построена в работах

[2-4] как регуляризованный вид системы уравнений Навье-Стокса (НС).

КГД система может рассматриваться как система уравнений НС, усреднен-

ная на малом пространственно-временном интервале, что приводит к сгла-

живанию (регуляризации) исходной системы уравнений.

КГД система уравнений в декартовых координатах в отсутствие внеш-

них сил и источников тепла может быть представлена в виде [3]:







(1)

i mj

t



u



(

) 





(2)

t





(

H)



(

)

(3)

t

Здесь  - плотность газа,uⁱ - компоненты его макроскопической скорости,

p - давление. Полная энергия единицы объема E и полная удельная эн-

тальпия H идеального политропного газа с показателем адиабаты  вы-

числяются по формулам

E u

/2 p /( 1),

H (Ep)/.

(4)

Вектор плотности потока массы ^imj определяется в виде



(u

w



(

u



(5)



Выражения для тензора вязких напряжений^ij и теплового потокаqⁱ

записываются как

k j





 u

(u







p) 



p  p

(6)



j i



(



u 



)

 



(7)

q

 u

(u



 pu



)

 

T .

(8)



120

И.А. Широков, Т.Г. Елизарова

Здесь   p/(( 1)) - внутренняя энергия единицы массы газа,^ijNS иq^iNS

- тензор вязких напряжений и тепловой поток в системе НС;  ,  и  - ко-

эффициенты сдвиговой и объемной вязкости и теплопроводности соответ-

ственно, T - температура газа.

Коэффициент сдвиговой вязкости  определим через температурную

зависимость [3]:

 

(T /T

)^

(9)

где



- вязкость газа при температуре

, 0 1 - показатель межмоле-

кулярного взаимодействия. Коэффициент объемной вязкости может быть

вычислен с использованием аппроксимационной формулы [3]

 ((5/3) ),

(10)

коэффициент теплопроводности вычисляется как

 /(Pr( 1)).

(11)

Коэффициент  , определяющий дополнительную диссипацию в КГД

алгоритме, для вязкого политропного газа имеет порядок характерного вре-

мени между столкновениями частиц газа [2-4]. Его величина связывается с

коэффициентом сдвиговой вязкости и может быть вычислена в виде [3, 4]

 /(pSc),

(12)

где Sc - число Шмидта, которое для газов близко к единице.

Для обеспечения устойчивости КГД алгоритма при моделировании

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

сипативного коэффициента  (12) путем включения в него слагаемого, за-

висящего от шага пространственной сетки и параметров течения в виде

 /(pSc)h /c,

(13)

где h - характерный размер пространственной ячейки, c - локальная ско-

рость звука,  - настроечный параметр, который в большинстве случаев

полагается постоянным числом порядка 1 [2-5].

Однако при больших сверхзвуковых скоростях течений введенной  -

вязкости вновь оказывается недостаточно. Как показывает практика чис-

ленных расчетов, в этих случаях для стабилизации численного алгоритма

возможно использовать наличие объемной вязкости, вводя искусственную

добавку в коэффициент (10) в виде

 ((5/3) )(h /c)p.

(14)

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

121

Величина регуляризирующей добавки здесь также определяется ло-

кальными параметрами и зависит от настроечного коэффициента . Такое

введение искусственной диссипации в КГД уравнения впервые использова-

но в [6, 7].

КГД алгоритм построен на основе системы (1)-(8) с диссипативными

коэффициентами (9)-(12) и искусственными добавками (13) и (14).

4. Численный алгоритм

Схема расчетной области изображена на рис.2. Набегающий воздуш-

ный поток направлен в сторону положительных значений оси x , вектор

скорости набегающего потока лежит в плоскости z  0 .

Нерегулярная тетраэдральная расчетная сетка строится посредством

свободно распространяемой библиотеки TetGen [8]. На рис.3 приведен об-

щий вид сетки, на рис.4 - фрагмент сетки в окрестности модели в сечении

z

. На поверхности модели элементы сетки представляют собой прямо-

угольные равнобедренные треугольники. Основная сетка имеет следующие

параметры: общее число узлов 298403, число тетраэдральных элементов

1574869, число узлов на поверхности модели 81312. Для исследования

влияния размера сетки на точность моделирования использовалась более

грубая сетка с такими параметрами: общее число узлов 137991, число тет-

раэдральных элементов 732815, число узлов на поверхности модели 36269.

Таким образом, дополнительная сетка приблизительно в 2 раза грубее, чем

основная. Как показали расчеты, не только число узлов, но и качество про-

странственной сетки существенно влияет на точность описания полей тече-

ния.

Для численного моделирования применялся существенно доработан-

ный программный комплекс, созданный А.А Свердлиным и Э.М. Кононо-

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

намических течений для тел произвольной формы с использованием тетра-

эдральных неструктурированных пространственных сеток.

Газодинамические параметры (плотность, скорость, давление, темпера-

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

раметров выбраны характерная длина (1м), плотность



и скорость звука

в набегающем потоке. Значения газодинамических параметров опреде-

ляются в узлах сетки. Конечно-разностная аппроксимация макроскопиче-

ских КГД уравнений строится методом контрольных объемов. Решение на-

чально-краевой задачи для сеточных аналогов КГД уравнений находится по

явной по времени конечно-разностной схеме. Пространственные производ-

122

И.А. Широков, Т.Г. Елизарова

ные аппроксимируются со вторым порядком точности, производные по

времени  с первым порядком точности.

Рис.2. Общий вид расчетной области.

Рис.3. Сетка в сечении z  0 .

Рис.4. Фрагмент сетки в окрестности модели в сечении z  0.

В КГД алгоритм включена искусственная диссипация с настроечными

коэффициентами  и  (13),(14). Коэффициент   1 считался постоянным,

а при вычислении коэффициента  использовалась следующая зависимость

от локального числа Маха: при Ma  2   5, при Ma  1.5   2 , в проме-

жутке между значениями числа Маха 1.5 и 2 коэффициент  линейно воз-

растает с числом Маха [6].

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

123

В начальный момент на входной границе задаются параметры набе-

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

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

На входной границе значения набегающего потока поддерживаются

постоянными. На выходной границе ставятся условия сноса, позволяющие

газу свободно покидать область. На твердой границе обтекаемого тела ста-

вятся условия прилипания (вектор скорости равен нулю), при этом исполь-

зуется дополнительное граничное условие КГД алгоритма: нормальные

производные давления и плотности равны нулю [3-5].

Как показали аналитические оценки и практика численных расчетов

[2-4], КГД алгоритм обладает устойчивостью при выполнении условия Ку-

ранта. Поэтому шаг по времени вычисляется следующим образом:_t 



/ c, где





0.1 - число Куранта, h - характерный локальный размер

пространственной сетки, c  T

- локальная скорость звука.

Расчеты проводились на многопроцессорном комплексе К-100 [10].

Программный комплекс обладает хорошей масштабируемостью и доста-

точно высокой эффективностью распараллеливания. В представленных

расчётах использовано 128 процессорных ядер.

5. Результаты моделирования

На рис.5-8 в сечении z  0 приведены контуры плотности для четырех

комбинаций числа Маха и угла атаки (Angle of attack, AoA), эксперимен-

тально и численно исследованных в [1]. Также изображены линии тока. От-

метим, что используемая реализация КГД алгоритма позволяет единообраз-

но проводить моделирование в области, целиком окружающей модель, в

том числе за хвостовой частью.

Рис.5. Ma  3, AoA  10^ .

Рис.6. Ma  3, AoA  20^ .

На рис.9-11 показано распределение давления по поверхности модели,

отнесенного к давлению в невозмущенном потоке, для варианта Ma 3,

AoA10^ . Рис.10,11 показывают проекцию поверхности модели на плос-

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

125

Сплошная линия соответствует расчету на основной сетке с общим

числом узлов 298403. Пунктирная линия изображает результаты, получен-

ные на дополнительной, более грубой сетке с числом узлов 137991. Симво-

лы показывают результаты, полученные при экспериментальном изучении

обтекания модели в [1]. Экспериментальные данные являются значениями

давления, полученными посредством датчиков, расположенных на поверх-

ности модели.

Треугольные символы соответствуют наветренной стороне модели

( y  0), квадратные - подветренной стороне ( y  0). В целом можно видеть

очень хорошее соответствие результатов, полученных в настоящей работе,

экспериментальным данным. Также видно, что результаты, полученные на

подробной сетке (сплошная линия), в целом ближе к экспериментальным

данным, нежели полученные на более грубой сетке (пунктир).

Рис.12. Ma  3, AoA  10^ .

Небольшое завышение максимального давления в носовой области мо-

дели по сравнению с экспериментальным обусловлено влиянием сетки, по-

скольку на более подробной сетке максимальное давление почти соответст-

вует экспериментальному. На наветренной стороне носовой части модели

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

значениями. Расхождения присутствуют в районе хвостовой части. Отме-

тим, что в экспериментальном исследовании обтекания модель закреплена

на кронштейне, который может влиять на картину обтекания в районе хво-

стовой части. Что касается подветренной стороны модели, то завышение

результатов моделирования по сравнению с экспериментальными данными

наблюдается в районе носовой части. В этой части присутствует вихревое

движение потока.

Отметим, что при анализе результатов моделирования наблюдалась не-

стационарность потока, особенно в области вихрей. Возможности КГД ал-

горитма при моделировании нестационарных течений также показаны в [5].

126

И.А. Широков, Т.Г. Елизарова

Рис.13. Ma  4 , AoA  20^ .

На рис.13 изображены одномерные профили давления на поверхности

модели в сечении z  0 , отнесенного к давлению в невозмущенном потоке,

для варианта Ma  4, AoA  20^ . Результаты получены на основной сетке с

числом узлов 298403. Символы показывают результаты, полученные при

экспериментальном изучении обтекания модели [1] (обозначения те же, что

на рис.12). По сравнению с вариантом Ma  3, AoA  10^ (рис.12), расхож-

дение между результатами моделирования и экспериментом больше. Осо-

бенно это видно в завышении максимального значения давления на носовой

части, а также в отсутствии скачка уплотнения на конической части модели

(x  0.08). Вероятно, тут проявляется излишняя вязкость, введенная в КГД

алгоритм, а также существенная нерегулярность используемой относитель-

но грубой сетки.

Рис.14, аналогично рис.10, изображает проекцию подветренной сторо-

ны ( y  0) модели на плоскость y  0 и линии тока на поверхности, для ва-

рианта Ma  4 , AoA  20^ . Для сравнения на рис.15 изображена аналогич-

ная проекция, взятая из [1]. На ней показаны линии тока, представляющие

собой результат обработки экспериментальных данных при тех же парамет-

рах Ma  4 , AoA  20^ . Видно, что экспериментальная картина линий тока

(рис.15) в общих чертах воспроизводится при моделировании на основе

КГД алгоритма (рис.14). Особенно это заметно в расположении парных

вихрей перед конической частью модели.

Рис.16 изображает проекцию подветренной стороны ( y  0) модели на

плоскость y0 и линии тока на поверхности для варианта Ma  4, AoA

20^ . В этой части поверхности не возникает вихревое движение, поэтому

картина течения практически симметрична, как и на аналогичном рис.11.

На рис.17 показан пример распределения расчетной области по 128 процес-

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

127

сорам, используемым в одном из вариантов расчетов, в сечении z  0 . Об-

щее число шагов по явной разностной схеме для каждого из вариантов со-

ставляет порядка

210

, требуемое машинное время около 20 часов.

Рис.14. Ma  4 , AoA  20^ .

Рис.15. Ma  4 , AoA  20^ .

Рис.16. Ma  4 , AoA  20^ .

Рис.17. Распределение по процессорам.

Авторы решают задачу моделирования в приближении уравнений НС с

учетом вязких слагаемых, при этом на поверхности тела стоят условия при-

липания. Таким образом, возникает вопрос о разрешении пограничного

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

прояснения этого вопроса в представленных расчетах приведены рис.18 и

u вблизи поверхности

модели в сечении передней цилиндрической части x  0.03 м для варианта

Ma  4 , AoA  20^ . На рис.18 также изображена расчётная сетка. На рис.19

показан одномерный профиль скоростиu_x при x  0.03 м, y  0, на кото-

рый нанесены маркеры, соответствующие границам ячеек сетки.

128

И.А. Широков, Т.Г. Елизарова

Рис.18. Ma  4 , AoA  20^ .

Рис.19. Профиль скорости вблизи поверхности.

Видно, что в пограничном слое находятся порядка 7 ячеек сетки. Та-

ким образом, при моделировании происходит частичное разрешение погра-

ничного слоя. Отметим, что на точность описания пограничного слоя влия-

ет не только количество точек внутри него, но и качество сетки. Практика

численных расчетов показывает, что тетраэдральная сетка не является оп-

тимальной для описания пограничного слоя. Тем не менее, даже на сетке

невысокого качества, использованной в настоящей работе, КГД алгоритм

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

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

В [12] проводится численное моделирование сверхзвукового обтекания

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

турбулентности Спаларта-Аллмараса и SST Ментера и проводится сравни-

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

настоящей работе структура течения (в том числе вихри) и распределение

давления по поверхности модели в целом аналогичны результатам, полу-

ченным в [12]. Однако структура фронтов ударных волн, образующихся в

окрестности модели при сверхзвуковом обтекании, разрешается в [12] луч-

ше. Отметим, что в [12] применена существенно более качественная и под-

робная расчетная сетка, содержащая 5963967 шестигранных ячеек (для

сравнения, в настоящей работе сетка состоит из 1574869 тетраэдров). Прак-

тика численных расчетов показывает, что качество сетки существенно влия-

ет на точность численного моделирования.

Одна из версий КГД алгоритма включена в качестве дополнительного

вычислительного ядра в открытый программный комплекс OpenFOAM [11].

Авторы благодарны А.Е. Луцкому за ценное обсуждение данной тематики.

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

129

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

E.M. Houtman, W.J. Bannink, B.H. Timmerman. Experimental and Computational Study

of a Blunt Cylinder-Flare Model in High Supersonic Flow. - Delft: Delft University of

Technology, 1995, report LR-796 .

B.N. Chetverushkin. Kinetic Schemes and Quasi-Gas Dynamic System of Equations.

Barselona: CIMNE, 2008.

T.G. Elizarova. Quasi-Gas Dynamic Equations. Dordrecht: Springer, 2009, IBSN 978-3-

642-0029-5.

Ю.В. Шеретов. Регуляризованные уравнения гидродинамики.  Тверь: Тверской го-

сударственный университет, 2016, 222 с.;

Yu.V. Sheretov. Regularized Hydrodynamic Equations.  Tver: State University, 2016, 222p.

elibrary.ru/item.asp?id=30097584 (in Russian).

Т.Г. Елизарова, И А. Широков. Регуляризованные уравнения и примеры их использо-

вания при моделировании газодинамических течений.  М.:-МАКС Пресс, 2017, 136с.

https://elibrary.ru/item.asp?id=29352202

T.G. Elizarova, I.A. Shirokov. Regularized equations and examples of their use in the mod-

eling of gas-dynamic flows.

И.А. Широков, Т.Г. Елизарова. Компьютерное моделирование обтекания модели

сверхзвуковым потоком вязкого сжимаемого газа на основе квазигазодинамического

алгоритма // Физико-химическая кинетика в газовой динамике, 2017, т.18, вып.2.

http://chemphys.edu.ru/issues/2017-18-2/articles/721

I.A. Shirokov, T.G. Elizarova. Computer simulation of the supersonic flow of a viscous

compressible gas around a model body on the basis of the quasi-gas-dynamic algorithm //

Physical-Chemical Kinetics in Gas Dynamics, 2017, v.18, iss. 2. http://chemphys.edu.ru/is-

sues/ 2017-18-2/ articles/721

T.G. Elizarova, I.A. Shirokov. Artificial dissipation coefficients in regularized equations of

supersonic aerodynamics // Doklady Mathematics, 2018, v.98, №3, p.648-651.

TetGen: A quality tetrahedral mesh generator. http://tetgen.berlios.de/

Т.А. Кудряшова, С.В. Поляков, А.А. Свердлин. Расчет параметров течения газа вокруг

спускаемого аппарата // Математическое моделирование, 2008, т.20, №7, с.13-22;

T.A. Kudryashova, S.V. Polyakov, A.A. Sverdlin. Calculation of gas flow parameters around

reentry vehicle // Math. Models and Comput. Simul., 2009, v.1, №4, p.445-452 .

10. K-100 System, Keldysh Institute of Applied Mathematics RAS, Moscow; Available at

http://www.kiam.ru/MVS/resourses/k100.htm

11. M.V. Kraposhin, E.V. Smirnova, T.G. Elizarova, M.A. Istomina. Development of a new

OpenFOAM solver using regularized gas-dynamic equations // Computers and Fluids,

2018, v.166, p.163-175.

12. В.E. Борисов, А.А. Давыдов, А.Е. Луцкий, Я.В. Ханхасаева. Численное исследование об-

текания модели космического аппарата // Препринты ИПМ им. М.В. Келдыша, 2017,

№130, 19с. doi:10.20948/prepr-2017-130. http://library.keldysh.ru/ preprint.asp?id=2017-130;

V.E. Borisov, A.A. Davydov, A.E. Lutskii, Ia.V. Khankhasaeva. Chislennoe issledovanie ob-

tekaniia modeli kosmicheskogo apparata // Preprinty IPM im. M.V. Keldysha, 2017, №130,

19s.; doi:10.20948/ prepr-2017-130. http://library.keldysh.ru/preprint.asp?id=2017-130.

Поступила в редакцию 04.03.2019

После доработки 30.04.2019

Принята к публикации 01.07.2019