Научная статья на тему 'БЕССЕТОЧНЫЙ АЛГОРИТМ РАСЧЁТА СВЕРХЗВУКОВЫХ ТЕЧЕНИЙ ВЯЗКОГО ТЕПЛОПРОВОДНОГО ГАЗА'

БЕССЕТОЧНЫЙ АЛГОРИТМ РАСЧЁТА СВЕРХЗВУКОВЫХ ТЕЧЕНИЙ ВЯЗКОГО ТЕПЛОПРОВОДНОГО ГАЗА Текст научной статьи по специальности «Математика»

CC BY
10
3
i Надоели баннеры? Вы всегда можете отключить рекламу.
Журнал
Труды МАИ
ВАК
Область наук
Ключевые слова
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ / БЕССЕТОЧНЫЙ МЕТОД / УРАВНЕНИЯ НАВЬЕ-СТОКСА / СВЕРХЗВУКОВОЕ ОБТЕКАНИЕ ТЕЛ / ТЕПЛОВОЙ ПОТОК

Аннотация научной статьи по математике, автор научной работы — Способин Андрей Витальевич

Изложен бессеточный алгоритм численного решения системы уравнений Навье-Стокса для вязкого газа. Рассмотрено решение модельной задачи сверхзвукового обтекания сферы. Полученные результаты показали хорошее согласование с эталонными значениями параметров газа на обтекаемой поверхности и приближённо-аналитическими выражениями для конвективного теплового потока от газа к поверхности и положения головной ударной волны.

i Надоели баннеры? Вы всегда можете отключить рекламу.

Похожие темы научных работ по математике , автор научной работы — Способин Андрей Витальевич

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

MESHLESS ALGORITHM FOR CALCULATING SUPERSONIC VISCOUS GAS FLOWS

Тhe article presents in detail the algorithm for the numerical solution of the Navier-Stokes equations [12] by the meshless method [8, 10]. The described method is used for numerical simulation of blunt bodies flow-around by supersonic viscous gas flow [1, 2]. Cartesian-grid-based immersed boundary method presented in previously published works was successfully applied to simulate such flows in the two-dimensional planar and an axisymmetric formulation [3]. Numerical studies of the gas thermal impact change on the streamlined surface while the highly inertial particle motion against the incoming flow along the symmetry axis were performed by dint of it[4, 5, 6]. Particles motion along various trajectories required gas-dynamics problems solving in the 3D formulation, for which the Cartesian-grid based methods application required too much computer memory. When solving gas-dynamics problems by the meshless method, a finite sets of nodes is being selected in the space. Approximation by the least square method is applied for spatial derivatives computing in computational nodes. The said approach is being used for convective and viscous fluxes computing. The convection fluxes are being computed by the AUSMPW+ method in conjunction with the MUSCL scheme and van Abada limiter. The system of equations time integration is being performed by the explicit Runge-Kutta method of the third order of accuracy. The flow-around surface is being represented by the model of isothermal wall with the specified temperature. Conditions of gas adhesion, as well as pressure gradient equality to zero, which modeling also employs the least square method, is realized on the said model. The meshless method gives the opportunity to compute the gas flow in the areas with complex geometry, it is simpler in realization herewith compared to the finite volumes method, since it does not require generation of the high-quality computational grid. It allows setting anisotropic distribution of nodes in space, which is of vital importance for qualitative resolution of the boundary layer near the flow-around surface. Along with the 3D realization, adaptation of meshless method for computing flat and axisymmetric viscous gas flows was performed. The software implementation of the described method is realized in the C++ programming language and employs the OpenMP technology for computations parallelization. A computational experiment was performed on modeling the sphere flow-around by the supersonic airflow at the Mach number of M = 3 and Reynolds number of Re = 105 was conducted for the said method verification. The spatial shadow patterns of the flow, pressure field and Mach number are presented. Comparison of gas parameters in the boundary layer, obtained by the meshless method with the computational results of continuous flow combined with the boundary layer equations is presented. Correspondence of the gas computational parameters on the flow-around surface to the reference data and approximate-analytical expressions is demonstrated. The heat flow value in the critical point coincides with the one calculated with the Fay and Riddell formula, and the heat flow distribution curve along the surface is close to the approximate-analytical one. The next stage in the development of the meshless method is planned to solve numerically the unsteady multiphase flows problem, in particular, to simulate the solid particles motion in the shock layer.

Текст научной работы на тему «БЕССЕТОЧНЫЙ АЛГОРИТМ РАСЧЁТА СВЕРХЗВУКОВЫХ ТЕЧЕНИЙ ВЯЗКОГО ТЕПЛОПРОВОДНОГО ГАЗА»

Труды МАИ. 2021. № 121 Trudy MAI, 2021, no. 121

Научная статья

УДК 519.63, 532.511

DOI: 10.34759/trd-2021-121-09

БЕССЕТОЧНЫЙ АЛГОРИТМ РАСЧЁТА СВЕРХЗВУКОВЫХ ТЕЧЕНИЙ ВЯЗКОГО ТЕПЛОПРОВОДНОГО ГАЗА

Андрей Витальевич Способин

Московский авиационный институт (национальный исследовательский университет),

Москва, Россия

spise@inbox.ru

Аннотация. Изложен бессеточный алгоритм численного решения системы уравнений Навье-Стокса для вязкого газа. Рассмотрено решение модельной задачи сверхзвукового обтекания сферы. Полученные результаты показали хорошее согласование с эталонными значениями параметров газа на обтекаемой поверхности и приближённо-аналитическими выражениями для конвективного теплового потока от газа к поверхности и положения головной ударной волны.

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

Финансирование: работа выполнена в рамках государственного задания Минобрнауки России, номер темы FSFF-2020-0013

Для цитирования: Способин А.В. Бессеточный алгоритм расчёта сверхзвуковых течений вязкого теплопроводного газа // Труды МАИ. 2021. № 121. DOI: 10.34759/trd-2021-121-09

MESHLESS ALGORITHM FOR CALCULATING SUPERSONIC

VISCOUS GAS FLOWS

Andrey V. Sposobin

Moscow Aviation Institute (National Research University), MAI,

Moscow, Russia

spise@inbox.ru

Abstract. ^e article presents in detail the algorithm for the numerical solution of the Navier-Stokes equations [12] by the meshless method [8, 10]. The described method is used for numerical simulation of blunt bodies flow-around by supersonic viscous gas flow [1, 2]. Cartesian-grid-based immersed boundary method presented in previously published works was successfully applied to simulate such flows in the two-dimensional planar and an axisymmetric formulation [3]. Numerical studies of the gas thermal impact change on the streamlined surface while the highly inertial particle motion against the incoming flow along the symmetry axis were performed by dint of it[4, 5, 6]. Particles motion along various trajectories required gas-dynamics problems solving in the 3D formulation, for which the Cartesian-grid based methods application required too much computer memory.

When solving gas-dynamics problems by the meshless method, a finite sets of nodes is being selected in the space. Approximation by the least square method is applied for spatial derivatives computing in computational nodes. The said approach is being used for convective and viscous fluxes computing. The convection fluxes are being computed by the AUSMPW+ method in conjunction with the MUSCL scheme and van Abada limiter. The system of equations time integration is being performed by the explicit Runge-Kutta method of the third order of accuracy. The flow-around surface is being represented by the model of isothermal wall with the specified temperature. Conditions of gas adhesion, as well as pressure gradient equality to zero, which modeling also employs the least square method, is realized on the said model.

The meshless method gives the opportunity to compute the gas flow in the areas with complex geometry, it is simpler in realization herewith compared to the finite volumes method, since it does not require generation of the high-quality computational grid. It allows setting anisotropic distribution of nodes in space, which is of vital importance for qualitative resolution of the boundary layer near the flow-around surface.

Along with the 3D realization, adaptation of meshless method for computing flat and axisymmetric viscous gas flows was performed.

The software implementation of the described method is realized in the C++ programming language and employs the OpenMP technology for computations parallelization. A computational experiment was performed on modeling the sphere flow-around by the supersonic airflow at the Mach number of M = 3 and Reynolds number of Re = 105 was

conducted for the said method verification. The spatial shadow patterns of the flow, pressure field and Mach number are presented. Comparison of gas parameters in the boundary layer, obtained by the meshless method with the computational results of continuous flow combined with the boundary layer equations is presented. Correspondence of the gas computational parameters on the flow-around surface to the reference data and approximate-analytical expressions is demonstrated. The heat flow value in the critical point coincides with the one calculated with the Fay and Riddell formula, and the heat flow distribution curve along the surface is close to the approximate-analytical one.

The next stage in the development of the meshless method is planned to solve numerically the unsteady multiphase flows problem, in particular, to simulate the solid particles motion in the shock layer.

Keywords: numerical simulation, meshless method, Navier-Stokes equation, supersonic flows around bodies, heat flux

Funding: the research was carried out within the state assignment of Ministry of Science and Higher Education of the Russian Federation (theme No. FSFF-2020-0013) For citation: Sposobin A.V. Meshless algorithm for calculating supersonic viscous gas flows. Trudy MAI, 2021, no. 121. DOI: 10.34759/trd-2021-121-09

Введение

В статье изложен алгоритм численного моделирования обтекания затупленного тела сверхзвуковым вязким газовым потоком [1, 2] в двумерной и трёхмерной

постановке. Работа является очередным этапом в рамках серии [3-6] численных исследований эволюции ударного слоя вблизи обтекаемой преграды и изменения тепловых потоков от газа к поверхности вследствие движения крупнодисперсной частицы. В натурных экспериментах [7] наличие крупной частицы в ударном слое приводило к кратному усилению конвективного теплового потока к обтекаемой поверхности.

Построенная ранее математическая модель на основе метода погруженной границы на декартовых сетках [3] позволила изучить симметричные течения с движением частицы вдоль оси симметрии. Для расширения спектра исследований потребовался менее требовательный к вычислительным ресурсам алгоритм. Выбор был сделан в пользу набирающих в последнее время популярность бессеточных [8-11] методов решения уравнений газовой динамики. При относительной простоте реализации, не требующей построения сложных вычислительных решёток, бессеточные методы обладают достаточной применительно к рассматриваемой задаче точностью, и в то же время позволяют существенно сократить затраты памяти компьютера благодаря анизотропии расположения расчётных узлов.

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

Математическая модель течения вязкого газа

Система дифференциальных уравнений Навье-Стокса в консервативных переменных в совокупности с уравнением состояния идеального газа описывает течение вязкого газа в трёхмерном пространстве [12]

^ дР (q) | ас (q) | дИ (q)_ дГ (q) < дСу (д) < дН (д)

0

тхх тху

ТХ2

тххи + тxyv + - Ях

дг 1 1 дх 1 ду д2 1 дх ду 1 д2

Г ри Л г р л г рw Л

ри 2 ри + р рuv рuw

д = рv , Р = рuv , с = 2 р\> + р , и = рvw , Fv =

рw рuw рvw рм!1 + р

к риН у V рvн , V р^Н у

с =

е =

0

'ух

'уу Туг

кТухи + ТууУ + Ту^ - Я у у

Р +1 (и2 + V2 + ), Н = е + р, р = рЯТ

Н =

о

'2х

Т2хи + TZyV + - Я 2

р(У-1) 2 ' Р

где г - время, р - плотность, р - давление, Т - температура, и, V, w - компоненты скорости газа, у - показатель адиабаты, Я - газовая постоянная. Компоненты тензора вязких напряжений:

тхх ^ М

^ди дv дw дх ду д2

, туу

2

3

^'дv ди дw Л ду дх д2

^д^ ди дv д2 дх ду

тху тух м

^ д и ду Л ду дх

, т

хх * хх

м

ди дм дх дх

тух тху м

^ ду дмл дх ду

; дТ дТ „дТ

Компоненты плотности теплового потока: = —А—, ду = -А—, = —А

дх у ду

дх

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

3

Л ГТ1 ГТ1*

м = м

Т V Т у

Т ~ , с Не

2 с где Т*= 273.15К, С* = 110,4К, м* = 0.0000178-

Т + С

.2

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

вязкости: А = ——, где С р - удельная теплоёмкость газа при постоянном давлении, Рг р

Рг = 0.72 - число Прандтля.

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

Численное решение системы уравнений в частных производных осуществляется на дискретном множестве расположенных в расчётной области точек (узлов).

Рис. 1. Распределение узлов в Рис. 2. Облака соседних вычислительных

центральном сечении расчётной области узлов

Аппроксимация заданного в пространстве поля скалярной величины ф=ф(х,у, х) выполняется согласно методу наименьших квадратов [10-11, 13]. Для каждого узла / и облака окружающих его соседних узлов у е С/:

л дф фу =ф +Ах1]~

+ Ау

дф

У

ду

Л дф

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

+ Ахи —

и дх

+ 0 (к2 ), Юу , ^ + Ау2 + Ах2

/ аи

= ху — х1, Дуц = у у — у%, АхЦ = х / — х/, ДФЦ =Фу — ф/,

гу

у

I

]еС/

ю

У

Дфу — Аху

дф дх

Ду

дф

у

ду

дф Ах/ } —

дх

\2

Ш1П

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

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

дф дх

л дф

= X а]АФ], "Г"

¡С ] ] ду

г ^

X Аф], ^ = X У АФ]-

а У]

X Йк Ахгк X Йк Ахгк АУгк X Йк Ахгк Аг1к

к^Сг к^Сг

X Югк Ахгк АУгк X Йк АУгк X Йгк АУгк Аггк

к^Сг к^Сг к^С^

X ®гк Ахгк Аггк X Йгк АУгк Аг1к X Йгк АгИк

к^Сг к^Сг к^Сг

АУг]

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

Используя полученные приближения, запишем уравнения газовой динамики в полу-дискретной форме:

дг

2 X[аrJ (% - % ) + Р (С] - С г ) + У ( Щ - Нг)'

^[а] (% - % ур] (с ] - о; уГ] (

н] - н

здесь %], С], Нг] - значения векторов конвективных потоков в середине отрезка, соединяющего узлы г и ], % = %(дг-), С = С(^), Н = н(^), = (дг),

с;=с; (я,.), н;=н; (я,.), =(qJ), с;;.=с; ), н;=н; ).

Вычисление невязких потоков

В ходе работы над программой решения системы уравнений Навье-Стокса бессеточным методом был реализован ряд известных модификаций схем расчёта конвективных потоков семейств НКЬ и AUSM. Наилучшее соответствие результатов численного расчёта эталонным данным было получено посредством схемы AUSMPW+ [14], описание которой и приводится ниже.

Реализация вычисления потоков , С у, Ну методом AUSMPW+ потребовала

двух проходов по набору точек. На первом этапе выполняется реконструкция векторов

консервативных переменных д = , = ¥+ посредством применения к

j

Т

компонентам вектора ¥ = схемы MUSCL с ограничителем van

Albada [11]:

Vij = Vi +

4

(1 - ksi)Aj +(1 + ksi)(vj - Vi)|, Aj = 2Arj • V щ - (vj- щ)

Vij = Vj K1'" ksJ )Aj +(l + ksJ )(vj - Vi)], A+ = 2Arj ^Vj- Vi) =

si = max

2Aj (v j -Vi ) + £

0

a-2 +(vj -vi) +£

sj = max

2A+ (vj -Vi ) + £ ' 2 A+2 +(Vj -Vi) +£

rj=

f \

Xj Xj yj - уг

К ZJ - 2г.

, ^Vn =

Г dy \ f

dx n

dy

dy n

dy

dz n V

X amn (ym yn )

X ßmn (V -Vn )

m<=C„

X Ymn (ym yn )

m^C„

1 -13 значения параметров k = —, s = 10

Для каждой пары узлов г и j расчётной области вычисляются минимальное

(min ( P-m, Pkm )) и коэффЩиент

значение давления Pmin = min

с i + w

fp = min ke{/,j}, mgC-

min

pkm pkm

- ' + V v pkm pkm у у

на всех интерфейсах между узлами г и j и их

+ - + непосредственными соседями, здесь p , p - величина давления, для q и q

соответственно.

+

На втором этапе с использованием значений рт{п, /р и векторов состояния я

и я- для пары узлов г и ] производится расчёт векторов конвективных потоков. Для обобщения представленного метода введём вектора п и N, зависящие от направления оси вектора потока:

N = (0 Лх Пу iz 0)T, П =

(1 0 0 )T для Fj (0 1 0 )T для G j ,

(0 0 1)T для H

jj

тогда нормальная проекция скорости 3 = V • п. Ниже подробно описан порядок расчёта

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

вектора конвективного потока К, по оси х, в этом случае 3 = и. Алгоритмы

У

вычисления векторов С ¡у, Ну аналогичны приведённому с учётом вектора п и

проекций скоростей, при этом принимаются 3 = V и 3 = w соответственно. В зависимости от направления вектора гу выбираются д^ и д^ :

q L

q+, Xj < Xj

, Xj ^ x j

q r

q

q , xj < xj

q

^ , Xj ^ X j

Вычисляются энтальпия, скорость звука и числа Маха:

H

normal

'hl- +Hr- м:Л

, cs

пГ-1 И 12 H normal

y +1

c 1/ = /2

С

max 3 cs )

,2

c

s

тах

(3, Cs )

, 3l + 3r > 0 , 3l + 3r < 0

3L д. 3R

ML =—, MR = —

С1/

CÏ2

Р

здесь Н = е л---полная энтальпия, vт = V - 3п - касательная компонента скорости.

Р

Рассчитываются коэффициенты

M+ =

( Ml +1)2

4

Ml + M

L

2

\Ml\< 1 \Ml\> 1

, M- =

_(Mr-1)2 ; MR, 1 . \MR\> 1

my2 = M+ + M-,

<

<

<

<

<

р

а=-

16

1 (Ыь +1)2 (2 -Ыь) + аЫь (ы1 -1) , №ь\< 1

1 (1 + sign (М^ ))

Р

Я

а=-

16

Рs = РЬ Рь + РЯ РЯ,

1 (1 -sign (Мя ))

0

г л ^ РЬ

1

Рs

Ш1П

s у

1,

Р Ш1П

Л2

V

Ш1П ( Рь, РЯ ),

0

г \ с РЯ

й=1 - /з.

-1

V Рs у

3

Ш1П

1,

Р ш1

\2

Ш1П

V

Ш1П ( Рь, РЯ )

, Р8 = 0

, Рs ^0'

, Р8 = 0

, Рs ^0

М

+

, К1 >1

1 (Мя -1)2 (2 + Мя)-аМя (мЯ -1)2 , \Мя\< 1

М- =

Мь + Мя-[(1 -й)(1 + /я )-/ь ] , ту > 0

М+й(1 + /ь) , ту^ <0

М--я-(1 + /я) , ту^ >0

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Мя+ М+-[(1 -я)(1 + /ь)-/я] , ту <0

|МЛ| >1

Выражение для расчёта конвективного потока % имеет вид

% = Мь ' с^ •ф ь + МЯ ' ^^ •ф Я + Рs •

N

<

3

<

3

<

<

<

<

т

Ф = (р ри ру рИ)

Расчёт вязких потоков

Для определения значений компонентов тензора вязких напряжений т и плотности теплового потока q в каждом расчётном узле I вычисляются частные производные по направлению температуры и компонентов скорости:

ди дх

X "у (иу-иг ), ду =\Рг} (иу-иг), ■ - , §. = X Угу (ту-Тг) .

1 ^

, 2 ( ~ ди ду д^ \ ( д и д V \

II 3 1 2— дх V г дУ /' дz г) , ■ . . , тху г =« кдУ ч-- / дх г У

> ■ •• , qz\г =-Л

дТ

дz

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

т^У ттУ

потоков р , Чг , Иг согласно приведённым выше соотношениям.

Интегрирование системы уравнений

Интегрирование системы дифференциальных уравнений по времени выполняется явным методом Рунге-Кутты третьего порядка:

дЧг

дг

Я (qг ) = (qг), Я (ч ) = 2 Х[агу (% %г ) + Ру (о у-Сг )+71] (Н уИ)

(чг )= X (%у-)+Ру (оу-оу Угг, (ну - НУ)

С (л\\

„« =яП -А к (ЧП), ч( 2)=3 яП+1 я!»-1А а Г

(1)

г

V У

яП+1 = 1 яП + 2 яР)- 2 Аг к

' (2 )л Я;

V у

Шаг по времени при явном интегрировании определяется согласно критерию Куранта

аг = ш1п

сбь ты

ах

V

шах /

шах /

Щ

Щ = X (ау% + ] + ]у + ] а + р] + У]

^С V

1 V

Щ = шах

' 3Л

4 у2

3, Р7

V У

X

Щ(а2 +рр] +у]) _ _

Рг+Р]

Р] =Р

ЛТ + Т] >

V 2 У

СБЬ = 0.5, Ут = 0.4

по Яое-усреднению векторов физических переменных для я+ и я

-

_ в + в +7Р+

Р = у[р Р , в = —-н^—, с ■

л/Р +4р+ \

М2 + V2 + н2

у(У- 1),

Л М Л/Р + ^ + Л/Р+ ~ V л/Р + ;+л!р+ Л н л/Р + Н ' \/Р М =--¡=^=-, V =--¡=^=—, н = * *

\l~p~+>/

Р

Р

Р

Моделирование плоских и осесимметричных течений

Представленный алгоритм несложным образом адаптируется для моделирования плоского течения в двумерной постановке, в этом случае параметры газа по оси £

дф

принимаются тождественно равными нулю н = 0, у] = 0, ——

г дг

= 0, а метод

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

т1

У

Д

у

Е Щк Ахгк Е Щк Ах1к АУ1к

к^С^ к^С^

Е Щк Ах1к АУ1к Е Щк АУ1к

к^С^ к^С^

Щу Аху

Щу АУу

Решение осесимметричной задачи может быть выполнено в цилиндрической системе координат, при этом система уравнений Навье-Стокса преобразуется к виду [12]:

дд2Б , дГ2б Б ) , дС2Б Б ) , дГ2Б (Я2Б ) , дС2Б (Я2Б ) , 1

--\---\---\— Кс =--1---1—

д? дх ду г дх ду г

' ри 1 ' ру 1 ' ру 1

д2 б = ри II С) 2 ри + р , С 2 Б = риу 2 > КС = риу 2

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

ру риу ру + р ру2

Р у \ риН у \ руН у \руН у

Г2Б =

0

'хх

сху

Тххи + тхуу - Цх

су

2 Б

0

'ух туу

\тхуи + тууу - Цу

К у

С

0

тху туу - твв \тхуи + тууу - Цу у

тхх ^ №

^ д и ду V дх ду г

*уу =-№

2 3

^^ду ди уЛ ду дх г

^ V д и ду г дх ду

здесь х - ось симметрии, ось у ей ортогональна, г - расстояние до оси симметрии.

Вычислительный эксперимент

Программная реализация представленных методов позволяет решать задачу моделирования обтекания тела сверхзвуковым вязким потоком газа в трёхмерной, а также двумерной плоской и осесимметричной постановках. Выполнена на языке программирования C++ для ОС Linux и Windows, использует технологию распараллеливания вычисления OpenMP [15, 16].

Для верификации представленных алгоритмов и их программной реализации решалась модельная задача сверхзвукового обтекания сферы радиусом Rs = 55 мм воздушным потоком. На входной границе расчётной области задаётся граничное условие первого рода, определяющее невозмущённый сверхзвуковой поток температурой T = 300K при числе Маха M= 3. Число Рейнольдса набегающего

потока Re= Ds = 105, здесь Ds = 2RS - диаметр сферы. На выходной границе М

области расчёта определены однородные условия Неймана.

Поверхность обтекаемого тела рассматривается как твёрдая изотермическая стенка с температурой Tw = 500 K, условиями непротекания и прилипания: v = 0,

T = Tw, —— = 0. Реализация граничных условий Неймана также основана на —n

аппроксимации частных производных — методом наименьших квадратов [17]:

—n

—т

—n

п-У(р = пх X a,i{(pj-(p,) + "y X ßu{(Pj-(Pi) + nz X Yi.i{(P.i-(Pi)

j^C, jeCi jeCi

Л/у = а1упх + Дупу +71упх ,

др

дп

X X щ(р]-(р! X щ

X -

др

р =

7 е С,'

дп

X %

где пх, Пу, п2 - компоненты вектора внешней нормали п в узле / на границе

поверхности, С; - множество его соседних узлов, не принадлежащих границе (см. рис. 3). Отсюда:

Т/ = тм>, Р/

X щр., Ел

р =~рЬ', и/ = у = м>/ = 0, е = р

ят

р (г-1)

Рис. 3. Узлы на границе и вблизи поверхности преграды

Рис. 4. Расположение узлов на сферической поверхности

Распределение вычислительных узлов во внешней области расчёта представляет собой прямоугольную декартову сетку с равномерным шагом. Вблизи окружающей обтекаемое тело области пространства узлы располагаются на поверхностях концентрических сфер, расстояние между которыми уменьшается по мере приближения к поверхности преграды (см. рис. 1). Минимальное расстояние между узлами задаётся внутри пограничного слоя, что необходимо для учёта вязких эффектов, так в критической точке на толщину пограничного слоя приходится не менее 50 слоёв вычислительных узлов. Узлы на каждой из концентрических сфер распределяются относительно равномерно согласно алгоритму, описанному в предыдущей работе по моделированию невязких течений [18, 19] (см. рис. 4).

В зоне стыковки сеток выполняется корректировка положения узлов внешней области при слабой обусловленности матрицы аппроксимации частных производных, вызванной слишком близким расположением соседних узлов. Схематическое расположение расчётных узлов в сечении области плоскостью 0Ху, проходящей через

центр обтекаемой сферы приведено на рис. 1.

На рис. 5 представлены полученные в вычислительных экспериментах теневые Шлирен-картины течения газа [20]. Штриховой кривой отмечено положение головной ударной волны, которое определяется приближённо-аналитическим выражением [21]:

X

(у ) = - Я, -А + Яс (МI -1)

1

1+

у

Я2 (М* -1)

-1

Яс = 1.143 • Я, • ехр

г л

0.54

(Мю-1)

1.2

, А = 0.143 • Я, • ехр

с \ 3.24

V М* J

Можно отметить хорошее соответствие между аналитическим и полученным численным методом пространственным положением головного скачка уплотнения. Некоторая зашумлённость теневой картины в целом характерна для методов семейства ЛиБЫ, а их использование объясняется меньшей погрешностью расчёта величины

конвективного теплового потока от газа поверхности по сравнению с методом ИЬЬ.

Рис. 5. Теневая картина течения (Шлирен)

30 60

а

Рис. 6. Коэффициент давления на поверхности сферы

На рис 6 приведено сравнение коэффициента давления ср на поверхности сферы,

полученного бессеточным методом, с эталонными значениями [22], наблюдается очень

хорошее соответствие результатов численного расчёта, коэффициент давления

вычисляется согласно формуле

с = 2 ср 2

Р - Рю

Рю IV

оо т оо

уМ{

ею V

Рю

- 1

где рю, рю, Vю, Мю - давление, плотность, скорость и число Маха набегающего

потока соответственно.

Поля давления газа и числа Маха в ударном слое представлены на рис. 7-8.

Рис. 7. Поле давления в ударном слое Рис. 8. Поле числа Маха в ударном слое

СПП _I_i_I_i_i_i_I_I_I_I_I_i_i_i_I_I

О 0.002 0.004 0.006 0.008 0.01 0.012

0.2

0.6

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

0.8

0.4

О1-1-1-

О 30 60

а

Рис. 10. Тепловой поток от газа к поверхности сферы

\

X

\

meshless

analytical approximation (Fay-Riddell, Polezhaev)

30

60

90

Рис. 9. Температура газа в пограничном слое

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

На рис. 10 приведены кривые теплового потока от газа к поверхности сферы. В качестве эталонной используется пригодная для инженерных расчётов простая зависимость, приближённо описывающая распределение теплового потока по поверхности сферы при высоких числах Маха в зависимости от угла а [24]:

0(а) = (0.55 + 0.45cos2а).

В критической точке эталонное значение конвективного теплового потока к поверхности тела с радиусом затупления Я8 вычисляется по формуле Фэя и Ридделла [25]:

с 1 л

/

дрк = а763|^^ ^ 06(Р0М0)04(Рм>Рм>)°Л(Н0 )

0.4,

у0.1

бу

|2

Ч^у0

гр _ гр

' Т° " Тда + 2С

Р

Р 0 = Р.

у

/ ? л-7 Г

(У + 1)М2

1

да

7-1

у + 1

V

2 уЫда-7 +1

7-1 Р0 Р^

' Р0 — ^Г , — Р0 , А _

ЛТ

0

лт

здесь: / — 1 для плоских течений и / — 0 для осесимметричных, энтальпия Н — С^рТ

Ср - теплоемкость газа при постоянном давлении, нижние индексы 0 соответствуют

параметрам газа в точке торможения, ^ - на стенке, да - набегающему невозмущенному потоку. Градиент скорости в критической точке рассчитывается по формуле:

1 2(Р0 - Рда)

I йУ) 0 А0

Отметим некоторые отклонения результатов расчёта теплового потока бессеточным методом решения уравнений Навье-Стокса от приближённо-

аналитической кривой, которые в проведённых расчётах не превышали 5%.

Заключение

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

Следующие ступени реализации проекта предполагают применение бессеточного метотода к расчёту движения частицы в потоке газа, а также реализацию приведённых алгоритмов для графических графических процессоров и вычислительных кластеров с целью расширения масштаба решаемых задач.

Список источников

1. Ревизников Д.Л., Сухарев Т.Ю. Гиперзвуковое обтекание затупленных тел в условиях атмосферы Земли и Марса. Сравнительный анализ математических моделей // Тепловые процессы в технике. 2018. Т. 10. № 1-2. С. 5-15.

2. Быков Л.В., Никитин П.В., Пашков О.А. Математическое моделирование процессов обтекания затупленного тела высокоскоростным потоком // Труды МАИ. 2014. № 78. URL: http://trudymai.ru/published.php?ID=53445

3. Винников В.В., Ревизников Д.Л., Способин А.В. Двухфазный ударный слой при обтекании тел сверхзвуковым запыленным потоком // Математическое моделирование. 2009. T. 21. № 12. C. 89-102.

4. Ревизников Д.Л., Способин А.В., Иванов И.Э. Изменение структуры течения под воздействием высокоинерционной частицы при обтекании тела сверхзвуковым гетерогенным потоком // Теплофизика высоких температур. 2018. Т. 56. № 6. С. 908913. DOI: 10.31857/S004036440003569-9

5. Reviznikov D.L., Sposobin A.V. and Ivanov I.E. Oscillatory flow regimes resulting from the shock layer - particle interaction // High Temperature, 2020, vol. 57, no. 2, p. 278283. D0I:10.1134/S0018151X20020169

6. Способин А.В., Ревизников Д.Л., Иванов И.Э., Крюков И.А. Колебания давления и теплового потока, индуцированные газодинамическим взаимодействием высокоинерционной частицы с ударным слоем // Известия высших учебных заведений. Авиационная техника. 2020. № 4. C. 108-115.

7. Holden M.S. et al. An Experimental Study of Particle-Induced Convective Heating Augmentation // AIAA Pap., 1976, no. 76, pp. 320. URL: https://doi.org/10.2514/6.1976-320

8. ^лстых А.И., Широбоков Д.А. Бессеточный метод на основе радиальных базисных функций // Журнал вычислительной математики и математической физики. 2005. Т. 45. № 8. С. 1498-1505.

9. Vasilyev A.N., Kolbin I.S., Reviznikov D.L. Meshfree computational algorithms based on normalized radial basis functions // International Symposium on Neural Networks, 2016, DOI:10.1007/978-3-319-40663-3 67

10. Sattarzadeh S., Jahangirian A. 3D implicit mesh-less method for compressible flow calculations // Scientia Iranica, 2012, vol. 19, no. 3, pp. 503-512. DOI:10.1016/j.scient.2012.04.013Y

11. Y. Wang, X. Cai, M. Zhang, X. Ma, D. Ren, J. Tan The study of the three-Dimensional meshless solver based on AUSM±up and MUSCL scheme // International Conference on Electromechanical Control Technology and Transportation, 2015. DOL1Q.2991/icectt-15.2015.52

12. Молчанов А.М. Математическое моделирование задач газодинамики и тепломассообмена. - М.: Изд-во МАИ, 2013. - 206 с.

13. Дринфельд Г.И. Интерполирование и способ наименьших квадратов. - Киев: Вища школа, 1984. - 103 с.

14. K.H. Kim, C. Kim, O.H. Rho. Methods for the accurate computations of hypersonic flows I. AUSMPW+ Scheme // Journal of Computational Physics, 2001, vol. 174, pp. 38-80. DQL10.1006/jcph.2001.6873

15. Антонов А.С. Параллельное программирование с использованием технологии OpenMP. - М.: Изд-во МГУ, 2009. - 79 с.

16. Малявко А.А. Параллельное программирование на основе технологий OpenMP, CUDA, OpenCL, MPI. - М.: Юрайт, 2021. - 135 с.

17. Hashemi M.Y., Jahangirian A. Implicit fully mesh-less method for compressible viscous flow calculations // Journal of Computational and Applied Mathematics, 2011, no. 235, pp. 4687-4700. DQI:10.1016/i.cam.2010.08.002

18. Способин А.В. Бессеточный алгоритм расчёта сверхзвуковых течений невязкого газа // Труды МАИ. 2021. № 119. URL: http://trudymai.ru/published.php?ID=159777. DOI: 10.34759/trd-2021-119-04

19. M. Deserno. How to Generate Equidistributed Points on the Surface of a Sphere, 2004. URL: https://www.cmu.edu/biolphys/deserno/pdf/sphere equi.pdf

20. Бодрышев В.В., Абашев В.М., Тарасенко О.С., Миролюбова Т.И. Интенсивность изображения, как количественная характеристика параметров газового потока // Труды МАИ. 2016. № 88. URL: http://trudymai.ru/published.php?ID=70428

21. Billig F.S. Shock-wave shapes around spherical-and cylindrical-nosed bodies // Journal of Spacecraft and Rockets, 1967, vol. 4, issue 6, pp. 822-823.

22. Любимов А.Н., Русанов В.В. Течение газа около тупых тел. - М.: Наука, 1970. ч. 1. - 287 с.

23. Ревизников Д.Л. Коэффициенты неизотермичности в задаче нестационарного сопряженного теплообмена на поверхности затупленных тел // Теплофизика высоких температур. 1995. Т. 33. № 2. С. 261-267.

24. Полежаев Ю.В., Юревич Ф.Б. Тепловая защита. - М.: Энергия, 1976. - 392 с.

25. Никитин П.В., Павлюк Е.А. Расчёт тепло-и массообмена на поверхности спускаемого космического аппарата // Труды МАИ. 2014. № 72. URL: http://trudymai.ru/published.php?ID=47266

References

1. Reviznikov D.L., Sukharev T.Yu. Teplovye protsessy v tekhnike, 2018, vol. 10, no. 1-2, pp. 5-15.

2. Bykov L.V., Nikitin P.V., Pashkov O.A. Trudy MAI, 2014, no. 78 URL: http://trudymai.ru/eng/published.php?ID=53445

3. Vinnikov V.V., Reviznikov D.L., Sposobin A.V. Matematicheskoe modelirovanie, 2009, vol. 21, no. 12, pp. 89-102.

4. Reviznikov D.L., Sposobin A.V., Ivanov I.E. Teplofizika vysokikh temperature, 2018, vol. 56, no. 6, pp. 908-913. DOI: 10.31857/S004036440003569-9

5. Reviznikov D.L., Sposobin A.V. and Ivanov I.E. Oscillatory flow regimes resulting from the shock layer - particle interaction, High Temperature, 2020, vol. 57, no. 2, p. 278-283. DOI:10.1134/S0018151X20020169

6. Sposobin A.V., Reviznikov D.L., Ivanov I.E., Kryukov I.A. Izvestiya vysshikh uchebnykh zavedenii. Aviatsionnaya tekhnika, 2020, no. 4, pp. 108-115.

7. Holden M.S. et al. An Experimental Study of Particle-Induced Convective Heating Augmentation, AIAA Pap., 1976, no. 76, pp. 320. URL: https://doi.org/10.251476.1976-320

8. Tolstykh A.I., Shirobokov D.A. Zhurnal vychislitel'noi matematiki i matematicheskoi fiziki. 2005, vol. 45, no. 8, pp. 1498-1505.

9. Vasilyev A.N., Kolbin I.S., Reviznikov D.L. Meshfree computational algorithms based on normalized radial basis functions, International Symposium on Neural Networks, 2016, D0I:10.1007/978-3-319-40663-3 67

10. Sattarzadeh S., Jahangirian A. 3D implicit mesh-less method for compressible flow calculations, Scientia Iranica, 2012, vol. 19, no. 3, pp. 503-512. D0I:10.1016/j.scient.2012.04.013Y

11. Y. Wang, X. Cai, M. Zhang, X. Ma, D. Ren, J. Tan The study of the three-Dimensional meshless solver based on AUSM±up and MUSCL scheme, International Conference on Electromechanical Control Technology and Transportation, 2015. D0I:10.2991/icectt-15.2015.52

12. Molchanov A.M. Matematicheskoe modelirovanie zadach gazodinamiki i teplomassoobmena (Mathematical Modeling of Problems of Gas Dynamics and Heat and Mass Transfer), Moscow, Izd-vo MAI, 2013, 206 p.

13. Drinfel'd G.I. Interpolirovanie i sposob naimen'shikh kvadratov (Interpolation and the method of least squares), Kiev, Vishcha shkola, 1984, 103 p.

14. K.H. Kim, C. Kim, O.H. Rho. Methods for the accurate computations of hypersonic flows I. AUSMPW+ Scheme, Journal of Computational Physics, 2001, vol. 174, pp. 38-80. DOI:10.1006/jcph.2001.6873

15. Antonov A.S. Parallel'noe programmirovanie s ispol'zovaniem tekhnologii OpenMP (Parallel programming using OpenMP technology), Moscow, Izd-vo MGU, 2009, 79 p.

16. Malyavko A.A. Parallel'noe programmirovanie na osnove tekhnologii OpenMP, CUDA, OpenCL, MPI (Parallel programming based on OpenMP, MPI, CUDA technology), Moscow, Yurait, 2021, 135 p.

17. Hashemi M.Y., Jahangirian A. Implicit fully mesh-less method for compressible viscous flow calculations, Journal of Computational and Applied Mathematics, 2011, no. 235, pp. 4687-4700. DOI:10.1016/j.cam.2010.08.002

18. Sposobin A.V. Trudy MAI, 2021, no. 119. URL: http://trudymai.ru/eng/published.php?ID=159777. DOI: 10.34759/trd-2021-119-04

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

19. M. Deserno. How to Generate Equidistributed Points on the Surface of a Sphere, 2004. URL: https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf

20. Bodryshev V.V., Abashev V.M., Tarasenko O.S., Mirolyubova T.I. Trudy MAI, 2016, no. 88. URL: http://trudymai.ru/eng/published.php?ID=70428

21. Billig F.S. Shock-wave shapes around spherical-and cylindrical-nosed bodies, Journal of Spacecraft and Rockets, 1967, vol. 4, issue 6, pp. 822-823.

22. Lyubimov A.N., Rusanov V.V. Techenie gaza okolo tupykh tel (Gas Flow Near Blunt Bodies), Moscow, Nauka, 1970, Part I, 287 p.

23. Reviznikov D.L. Teplofizika vysokikh temperatur, 1995, vol. 33, no. 2, pp. 261-267.

24. Polezhaev Yu.V., Yurevich F.B. Teplovaya zashchita (Thermal protection), Moscow, Energiya, 1976, 392 p.

25. Nikitin P.V., Pavlyuk E.A. Trudy MAI, 2014, no. 72. URL: http://trudymai.ru/eng/published.php?ID=47266

Статья поступила в редакцию 22.10.2021; одобрена после рецензирования 29.10.2021; принята к публикации 21.12.2021.

The article was submitted on 22.10.2021; approved after reviewing on 29.10.2021; accepted for publication on 21.12.2021.

i Надоели баннеры? Вы всегда можете отключить рекламу.