CC BY
13Russian Journal of Nonlinear Dynamics, 2022, vol. 18, no. 1, pp. 61-82. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd220105
NONLINEAR PHYSICS AND MECHANICS
MSC 2010: 35Q30 76-10 76N15
Validation of RANS Turbulence Models
for the Conjugate Heat Exchange Problem
A. A. Chernova
This paper addresses problems of mathematical modeling of heat exchange processes in the pre-nozzle volume of a solid propellant rocket engine with a charge with starlike cross-section and a recessed hinged nozzle. Methods of mathematical modeling are used to solve the quasi-stationary spatial conjugate problem of heat exchange. An analysis is made of the influence of RANS turbulence models on the flow structure in the flow channels of the engine and on the computed heat flow distributions over the surface of the recessed nozzle. Methods of mathematical modeling are used to solve the quasi-stationary spatial conjugate problem of heat exchange. Results of validation of RANS turbulence models are presented using well-known experimental data. A comparison of numerical and experimental distributions of the heat-transfer coefficient over the inlet surface of the recessed nozzle for the engine with a cylindrical channel charge is made for a primary choice of turbulence models providing a qualitative agreement between calculated and experimental data. By analyzing the results of numerical modeling of the conjugate problem of heat exchange in the combustion chamber of the solid propellant engine with a starlike channel, it is shown that the SST k — w turbulence model provides local heat-transfer coefficient distributions that are particularly close to the experimental data.
Keywords: solid propellant rocket engine, recessed nozzle, mathematical modeling, conjugate heat exchange problem, RANS turbulence models, heat-transfer coefficient
1. Introduction
There exists a large variety of approaches and models [8, 9, 12, 42, 44] for numerical description of turbulent flows. The application of each of them requires not only using appropriate forms of the system of conservation equations [8, 9, 12, 44], but also a different degree of resolution of the grid domain and, as a consequence, different computational resources.
Turbulence models falling into the class of Reynolds (both low-Reynolds and high-Reynolds) models [2, 13, 24, 27, 39], namely, the RANS and URANS turbulence models, are in most
Received September 24, 2021 Accepted February 03, 2022
Alena A. Chernova [email protected]
Kalashnikov Izhevsk State Technical University ul. Studencheskaya 7, Izhevsk, 426069 Russia
common use [44] since only minor requirements are placed on the discretization of the calculated area with respect to space and since they require low computational capacity. Application of such models for description of the turbulent flow involves the Reynolds and Favre averaging of the initial system of conservation equations, followed by application of a set of equations for taking into account turbulent viscosity, vorticity, and dissipation of turbulent energy. The family of k — e models may be regarded as the most frequently used [9, 44]. It includes the standard k — e model [2], the low-Reynolds k — e model [13]), the family of k — w models [2], the Spalart-Allmaras (SA) model [39] and various combinations of these models. The diversity of Reynolds models implies the application of a model that captures specific physical processes particularly well and takes into account not only viscosity and velocity, but also special features of the forming flows (prevalence of shear or near-wall flows, flows similar to jets, etc.). However, the use of a specific turbulence model (including a model from the RANS/URANS family) for a concrete physical problem requires preliminary testing.
Much attention is given in the literature to problems concerning the applicability of turbulence models for a wide range of problems. For example, in [29] an estimate is given of the influence of RANS turbulence models on the accuracy in modeling the process of the supersonic outflow of an ideal gas from a conical nozzle. No essential differences have been found in the structure of the supersonic flow of the gas expanding in the nozzle using the models considered. Ref. [40] deals with the problems of applying various turbulence models to simulate gas flows and to transfer admixtures in elements of urban areas. The problems of the applicability of RANS turbulence models with different degrees of elaboration for modeling the processes of internal ballistics in the combustion chamber of a gas turbine engine are examined in [12]. In particular, it is shown that the 1996 Menter model [24] (referred to as SST in what follows) captures more accurately the kinetics of the flow in the flow channel of the combustion chamber built according to the design under consideration. Aspects of the applicability of differential of turbulence models to the resolution of axisymmetric jets are examined in detail in [19], and the problems of modeling the turbulent transfer in high-temperature jets are treated in [18, 35]. But very little attention has been given in the recent literature to the problems of extending validation studies to similar physical objects built according to other designs.
Thus, the choice of an adequate turbulence model is determined by the type of flow, the specificity of the problem to be solved, and by the requirements placed on the accuracy of resolution of the phenomena of turbulent transfer. This paper addresses the problems of testing the applicability of RANS turbulence models to the modeling of intrachamber processes in the combustion chamber of a solid propellant rocket engine with a recessed nozzle.
2. Physical problem statement
Consider the operational processes proceeding in the combustion chamber of the solid pro-pellant rocket engine. When the exploder (pyrocartridge) comes into action, a hot gas wave propagates and the pressure increases instantaneously in the combustion chamber [17, 21]. The pyrocartridge itself is a volume source of heat release [1]. All this results in the solid propellant being ignited [21].
In the case described above, since the flow from the hot gas igniter is nonuniform, the spread of the flame front over the surface of the charge channel is characterized, on the one hand, by some nonuniformity and, on the other hand, by a small time interval. The process of solid propellant
ignition is generally regarded as instantaneous [22], and is not considered in the general form of studies of operational processes in the combustion chamber of the solid propellant rocket engine.
The initiation of the operational processes in the solid propellant rocket engine is accompanied by an additional increase in the pressure up to calculated values. The increase in the pressure leads to destruction (pop-out) of the plug, to the formation of a jet in the nozzle of the solid propellant rocket engine and thus to the rated operation conditions of the engine. The operation of the combustion chamber of the solid propellant rocket engine under the rated operation conditions is regarded, according to [22], as stationary from the point of view of changes in integral characteristics and as quasi-stationary from the point of view of the intrachamber processes. For this reason, the approach that involves fixation of instantaneous changes in the intrachamber parameters for some positions of the combustion vault is in common use [1, 3-6, 21, 22].
According to [22], intrachamber processes in the solid propellant rocket engine include the following interrealated physical processes:
• triggering of the pyrocartridge (combustion processes in the igniter);
• spread of combustion products over the combustion chamber of the solid propellant rocket engine;
• heat transfer from the combustion products to the propellant charge;
• heating-up and subsequent ingition of the solid propellant charge;
• nonstationary combustion of the charge;
• quasi-stationary combustion of the charge;
• FSI interaction of the actuation gas with the housing of the combustion chamber and the charge channel;
• change in the geometry of the flow channels of the combustion chamber of the solid propel-lant rocket engine.
Generalization of all physical processes proceeding in the combustion chamber of the solid propellant rocket engine within the framework of a unified mathematical model remains extremely difficult [1, 22] and in most cases unpractical. For this reason, it is common practice [22] to regard intrachamber processes as independent processes when it comes to an object or a time interval of the engine's operation. Specifically, the following classes of subproblems for description of all intrachamber processes are singled out:
• Ignition of the solid propellant rocket engine.
• Intrachamber processes under nonstationary and transient-state operation conditions of the solid propellant rocket engine, including nonstationary combustion of the charge.
• Internal ballistics of the solid propellant rocket engine, including quasi-stationary combustion of the solid propellant charge, when the engine operates in a stationary area.
• Stress-strain state of fuel elements of the charge channel and of the housing of the solid propellant rocket engine.
Thus, we will consider the flow of solid propellant combustion products in the flow channels of the combustion chamber of the solid propellant rocket engine. The description of gas-dynamical and thermophysical processes in the combustion chamber of the solid propellant rocket engine is based on conservation equations [16, 21, 22] averaged for the resolution of turbulent flows and velocity pulsations, and on the solution of the equation of nonstationary thermal conductivity for calculation of thermal conditions in the combustion chamber. The mass flux from the solid propellant combustion surface is nonuniform and nonstationary [17, 21], but can generally be represented by a uniform (average) consumption distributed over the entire area of the surface of the burning vault or can be defined as normal to the surface of the injection channel [42-44]. The combustion products are heat-conducting compressible gas and have an initial temperature equal to the combustion temperature of solid propellant.
3. Mathematical models
The conjugate problem of heat exchange in the pre-nozzle volume of the combustion chamber of the solid propellant rocket engine with a recessed hinged nozzle is considered. In [5, 6, 23, 37, 38, 42, 43] it is shown that a spatial turbulent (Re ^ 105) flow is obtained in the pre-nozzle volume. Since the gas velocity in the area of the nozzle entry (Fig. 1) reaches values corresponding to M ^ 0.7, [37, 38], and the gas velocity in the critical section of the nozzle, to the velocity of sound, the flow of combustion products is compressible.
Fig. 1. Generalized scheme of the solid propellant rocket engine with a recessed nozzle
The following assumptions are made in the mathematical description of intrachamber processes in the combustion chamber of the solid propellant rocket engine:
1. Assumption of the properties of solid propellant combustion products.
Since some special-purpose modifications of the solid propellant rocket engine, for example, solid propellant rocket propulsion systems, use fuel with a small content of condensed particles, the amount of the k-phase in the combustion products is assumed to be small and is not taken into account, and the gas mixture is considered to be homogenized and isotropic. It should be noted that, according to [7], the contribution of radiant heat transfer exchange in the flow with solid particles with volumetric heat release from the k-phase in the range of the recessed nozzle is 20-25 % in the subsonic region and 45-60 % in transand supersonic regions.
2. Assumptions of the pattern of solid propellant combustion.
The process of combustion of the charge of solid propellant is not taken into account and is replaced with distributed injection. The displacement of the combustion vault is not taken into account. The injection velocity from the mass flux surfaces is calculated as Pfuf — = Pg Ug.
3. Assumptions of the smoothness of the wall surface.
The microrelief of the wall surface of the regions under consideration is not taken into account.
4. Assumption of the structure of the walls of the combustion chamber.
The multilayer structure of the walls of the combustion chamber (Fig. 2) is considered. It consists of a heat-resistant coating layer (in the area of interaction with the flow of combustion products) and a metal layer (in the area of contact of the housing of the combustion chamber of the solid propellant rocket engine with the environment).
Indestructible "barrier" heat-resistant coatings consisting of refractory nonablative materials are considered.
Fig. 2. Scheme of the two-layer wall of the combustion chamber: the red line indicates the area of contact with the actuation gas (the boundary between the combustion products and the heat-resistant coating); the blue line denotes the external wall of the combustion chamber
5. Assumption of the pattern of heat exchange of the external surface of the last layer of the housing of the plant with the environment.
To describe the process of heat exchange between the external wall of the combustion chamber and the environment, use is made of the model of natural convection, which corresponds to the test conditions of power plants on a bench [38]. That is to say, the influence of external aerodynamics, including the aerodynamical convective heating of the housing, reemission of the walls etc., as well as the influence of changes in external conditions on the thermal conditions, is not taken into account.
The ambient temperature is taken to be constant and is defined in accordance with the air temperature during scavenging (300 K), and the heat-transfer coefficient a is determined from the dimensionless equation Nu — 0.15 • (GrPr)0-33.
Thus, the flow of combustion products in the flow channels of the combustion chamber of the solid propellant rocket engine can be described in terms of the model of viscous compressible
heat-conducting gas:
dU „ „ „ „ p— = pF-Vp + V-P, dt
di U2
P— ( c,„T +
dt
^ p = pRT
2
(3-1)
pFU + V • (PU) -V- q,
where p is the gas density, p is the pressure, U is the velocity vector, P is the stress tensor, F is the volume force, T is the temperature, q is the vector of the heat flow density, and R is the specific gas constant.
The thermal state of the walls of the combustion chamber of the rocket engine is determined from the equation of thermal conductivity:
pc^ = V(AV(T)),
(3-2)
where c is the thermal capacity and A is the coefficient of heat conductivity.
To construct the averaged forms of the conservation equations [31, 32], it is convenient to represent the system of equations (3.1), taking into account the equation of thermal conductivity (3.2), in the following divergent form:
{ dp dpUj dt dx„•
0,
dplTl dpUiUj
dt
dpE ~df
dx,
dp d
dx, dx, ^ j
p
+
dpEUj dpUj dUT
dx,
dx,
p = pRT,
c dT _ d xdT ^C dt dxi dxi '
rJU1 + dU1 dxj dxi
+ F.IJ.
dxi dxi
2 dU,
3" dxklj
Sij + F,,
(3-3)
where E = cvT + 0.5Uf is the total specific energy, t^ = 2pS^ — ^P^-^ij is the viscous stress
1 fdUi , dUj
tensor, S^ = ^ ^^ + -Q^rj is the tensor of deformation velocities, p is the dynamical coefficient of viscosity, and Sj is the Kronecker symbol.
We will perform the Favre averaging of the initial system of equations (3.3) [15], by transition to summation of the averaged and pulsation parameters, as \I> = \ir + \I>// (^ = p^/p, p^, p being
the Reynolds averaged parameters). Then the averaged system of equations takes the form d~p cfpLL
dt ^ dx-
0,
SrpLL , c7pUlU1 dp d ( , -
= -Ö--1- "ö [Ta + 7V,;,- + t ,;,
dt
+
dxj
"a?, ar,Vij tijJ
dXj dXj
cfpE cfpEU- _ dpUj _d
dt dXj
p = pRT, _ dT ,d2 T
~pcm- =
Ui (Ti:i + Ttij)] + ^ (<?, - qtj) + FjUj,
where r,, = a
Wi
9Xj
dUj
dx.
— |li^P-bij is the viscous stress tensor averaged in the sense of Favre,
Q]J dU \
Ttij = l-k I dx1 + 7ki~ J ~ Iis averaged turbulent stress tensor, ¿tt is the coefficient of
dynamical turbulent viscosity, k is the kinetic energy of turbulence, e is the turbulent dissipation, and qtj is the turbulent heat flow.
Since the number of unknowns increases, the averaged system of equations (3.4) loses closed-ness. Because of this, and to define a number of new variables, it is necesary to introduce additional equations defining the turbulent parameters, in particular, the turbulent viscosity. The methods of closing the averaged system of equations are determined by the turbulence model employed. Use is made of the following RANS turbulence models:
1. k — e. In this model, the equations of motion of the continuum (V2: environment) are supplemented with fluctuation terms determining the transfer of the kinetic energy of turbulence and the viscous dissipation of turbulence. When the k — e model is used, arbitrary small values are assigned to the turbulence characteristics on the permeable surface of the channel, and on the impermeable wall these characteristics are found by the method of near-wall functions. According to [20], the standard k — e turbulence model is defined by the following relations:
Bph + ^W) = 9
dt dpe ~dt
dx.
+
dx
dx.
d / rr ^ 9
(pUke) =
dx.
at = f,
cupk
at
ah J
at
ß + —
2
dk
dxk de dx.
+ Sk
+ Se,
k
e
2. k — w. This turbulence model is a correction of the model of D.Wilcox [46] taking the Kolmogorov hypothesis [14] into account, and is based on the solution to the equation of
the velocity of dissipation of the turbulent energy (w). This model has from the outset been aimed at the resolution of turbulent flows in boundary and shear layers and in flows
with low Reynolds numbers [44]. Problems of correction and adaptation of the model are
addressed in [9]. The model is defined by the equations
dk ^ ^ dk _ ^ dUi p*^ ^ dt 3 dXj 13 dXj dx
o
, * s dk
du du dt 3 dx
3
k 13 dXj dXj
k w '
(v + aut )
dw dx■
v
t
3. The SST model. This turbulence model, also known as BSL and as the Menter SST model [24-26], belongs to the class of two-band models, i.e., it uses equations of various turbulence models (in this case, k—e and k—w) [24] by incorporating the weight functions F1 and F2 into the formulation of the equations of the transfer of turbulence energy. Thus, the resulting model is described by the equations
dk ~dt
+ pUi^r- = pk ~ P*kuj +
' dx
3
JL
dx■
dw TT dw ~ n 2 d -rrr + Uj— = aS2 — ¡3u) + -—
dt 3 dx dxJ
dw vt =
. dk . . dw
^ . 1 dk dw
a1k
max(a1w, SF2) '
4. k—u SST. The Menter two-band turbulence model [11, 27], or the y—Re0 model, which is an evolutionary continuation of the SST turbulence model. In contrast to the initial Menter model [24], the equations of the updated k — u SST model include two more additional transfer equations and a set of empirical constants. The k — u SST turbulence model can be defined as
d(pl) d{pUf!)
dt
dx■
d
= P~fl ~ Ejl + P~f2 ~ E~(2 + g^r
p +
th
c?7 dx■
d ( pRee ~df~
+
d ( pUj Red
dx3
dRe$t
dx3
5. RNG k — e. The RNG-model of subgrid viscosity is an updating of the standard k — e model
based on renormalization theory [30]. Compared to the k — e model, the renormalized model
has been supplemented with e, another term in the equation which makes it possible to
enhance the accuracy of calculation of high-Reynolds flows. And incorporation of analytic
formulae of turbulent Prandtl numbers and effective viscosity into the model allows the
application of RNG k — e to be extended to low-Reynolds flows. The RNG k — e turbulence
model can be defined by the following equations:
d d d -{pk) + -{pklh) = —
1 J
d d d Ft(pe) + ^{peUÙ = —
1 J
C*2s = C2e +
a \ dk
'■'■<+-J to,
at +
de
ae ! dx■
eJ
1 + ßn3
Tj = S-
e
S = (2Sij Sij )1/2, at = Pcu — -
+ Pk - Pe,
ee + Cl eT^k — C-iePji
e
6. The Spalart-Allmaras model. This model can also be called the eddy viscosity transport equation (EVTE) and is one of the most widespread one-parameter turbulence models [36]. It is employed in modeling the external aerodynamics [39]. The transfer equation of the updated Spalart-Allmaras turbulence model has the form
d r d r ( T \ 2 1 d
^ + Cbl[1 ~ft2]s" ~CwlL u)
vi = rLi-
(v + r)
d r dx
k
+Cnfdv a
dx
2
All turbulence models presented above are supplemented with near-wall functions using the method due to Launder and Spalding [41]. In the area of the logarithmic layer the tangential velocity near the wall is related to the wall's shear stress, , by a logarithmic relation. The logarithmic dependence for the near-wall velocity is given as follows:
ut=l^ = \My+)+C1 (3.5)
Ut k
where
a vP
where u+ is the velocity near the wall, Ut is the known velocity that is tangent to the wall at distance Ay from the wall, y+ is the dimensionless distance from the wall, is the shear stress on the wall, k is the von Karman constant, and C is the logarithmic constant defined by roughness (for all the problems under consideration it is C = 0).
To solve the problem of singularity of the classical logarithmic dependence for the Launder-Spalding model of the near-wall velocity [41] (3.5), which arises as Ut ^ 0, for the k — e, RNG k — e and Spalart-Allmaras turbulence models, an alternative scale of velocities u* will be used instead of uT in the logarithmic region:
u* = C1/4k1/2. (3.7)
Then uT can be defined as
Ut
uT
ln(y8)
The absolute value of shear stress on the wall is
(3.8)
= pu uT
(3.9)
where y* is defined as y* = (pu*Ay)/p, and u* is defined as above. To rectify the discrepancy in the solution of the near-wall functions [28], use is made of the method of scalable near-wall functions or the method of restricting the values of y*, which is used in the logarithmic formulation, to the lower values — y* = max(y8, 11.06), where 11.06 is the value of y* at the point of intersection between the logarithmic and the linear near-wall profile. The calculated value of y* does not decrease below this limit. Thus, all points of the grid are located outside the viscous sublayer, and all small discrepancies of the grid are eliminated. The boundary condition for the velocity of dissipation e is given by the following relation valid in the logarithmic region:
£ =
pur
y* k
,3/4
(3.10)
For the SST, SST k — w and k — w turbulence models, an adapted formulation of near-wall function is applied which involves a successive refinement of the grid taking into account y+ from coarse grids which do not resolve the viscous sublayer to thin grids with the node lying inside the viscous sublayer. In this case, the equations read as
• The flow for the momentum equation (Fu):
Fu = —puT u*,
with
where
u
\
uT
AU
Ay
+ (y^fc)4, uT = y«)4 + «)4,
I fi P
AU
Ay
and uT = —
U
ln(y+) + C
The flow for the k-equation:
The expression for the w-equation:
Fh = 0.
=
u*
1 u*
a1 ky alkvy+
where ujs = ^^y (Ay being the distance between the first mesh point and the second mesh point. To achieve smooth mixing and to avoid cyclic convergence behavior, the following
formulation is chosen: = 1 + i Ay = y2 — y\-
4
k
It should be noted that the constants for all the turbulence models are defined in the standard way [44].
The discretization of the governing equations is performed by the finite-volume method, which is a grid method for solving nonlinear partial differential equations. The control volume (Fig. 3) is constructed around each mesh node (a tetrahedral grid element is considered as an example) using a median (defined by the lines connecting the centers of the ribs and the centers of the elements (facets) around the node).
Fig. 3. Construction of the control volume
The construction of the control volume method using, for example, the equations of conservation of the mass, momenta and some passive scalar can be represented as follows: 1. Finding the sought-for values at the center of control volume C:
J (r - rC) dV = 0.
Vn
2. The initial partial differential equations
dp d . TT. + —(pUj)
dt dx
0,
3
§№)+=-J|+ir u ( üj+if
dt
dx„•
dx„•
(3.11)
(3.12)
(3.13)
are integrated with respect to each control volume.
3. The application of the Gauss - Ostrogradsky theorem leads to the form
^ J pdV + J pUj diij = 0,
V S
(3.14)
d f f f f f d^U
pUi dV + / pUj Ui diij = — P diij + I ß
dt
j p diij I j /1 ( j diij + j SlT, dV, (3.15)
V
S
S
S
V
d 1 p(j)dV + f puj(f) diij = f r ( J-^- J diij + f S^dV,
dt
dx
(3.16)
V
S
S
3
V
where V and S are the volume and surface regions of integration.
4. Next, it is necessary to discretize volume and surface integrals. The volume integrals are discretized within the limits of each sector of an element (Fig. 4) and are summed up in the control volume to which the sector belongs. The surface integrals are discretized at the points of integration ipn (Fig. 4).
The point of integration
The center of the element
Fig. 4. Grid element
After discretization of the volume and surface integrals, the discretized equations take the
form
V
P-P At
+ Y1rh ip = 0
(3.17)
ip
V
PiUj-p°U? At
( ( dU- du3
+ ™ip(ui)ip = ¿^(p^nùip + 2^ [a h^1 +
ip
ip
ip
dx3 dxl 3 i
An, +SLTV,
ip
At
ipYip dx~ 3
ip op \ 3
+ V>
(3.18)
(3.19)
ip
where mip = (pUjAni)ip, V is the control volume, At is the time step, An, is the discrete external surface vector, the subscript ip denotes an estimate at the point of integration, the summation is over all points of integration of the control volume, and the superscript " 0" denotes the prededing time step.
When passing from the initial computational region to its discretized form, a replacement of the surface integrals is performed by summation over the facets of the control cells. In this case, the form of the approximation of the integrand is determined by the numerical scheme used.
The gradient in the control cell can be defined as
V
(3.20)
where Vi is the size of the cell i, sij is the facet between the cells i and j of the grid, nij is the unit vector of the external normal to the facet sij, and 0ij is the value of 0 on the facet ij.
The gradients on the facet between the cells are calculated by linear interpolation of the gradients from the cells i and j — V(0)j -.
The divergence is calculated from the formula
V
>,U, )
iv
s n • • l] l]
Uij 0
I] Y l] ■
The approximation of diffusion terms is constructed as follows:
1
s - -i] i]
V0.
l]1
where rij- is the diffusion coefficient on the facet ij.
To discretize the convective terms, the following scheme of central differences is applied:
0lp = 0up + ßV0
'up
(3.21)
where 0ip is the sought-for convective flow defined in terms of the values of the convective flow 0 at the nodes of the computational grid, 0up is the value of the convective flow in the node lying against the flow, r is the vector from the grid's angle lying against the flow to the point where the coefficient fi & 1 is determined locally with minimization of the local oscillations of solutions and is calculated separately for each component of the vector, and V0 is the value of the convective flow in the node lying against the flow.
The approximation of the diffusion terms of the conservation equations is constructed using shape functions. For example, for the derivative in the direction x at the point of integration ip:
00 dx
ip
ydNn dx
(3.22)
ip
The summation is taken over all shape functions for an element, and the derivatives of the shape functions with respect to Cartesian coordinates can be expressed in terms of their local derivatives using the Jacobi transformation matrix:
(3.23)
An estimate of the gradients of the shape function is constructed on the basis of trilinear interpolation. To resolve compressible flows, the product of the density and the advective velocity in the continuity equation is linearized:
~dN~ dx dy dz -1 'dN'
dx ds ds ds ds
dN dx dy dz dN
dy dt dt dt dt
dN dx dy dz dN
dz du du, du du,
(pU)nA & pnU0A + p0UnA — p°U0A, where the superscripts n and 0 denote the current and the preceding iteration.
(3.24)
r
The value of pn is linearized in the pressure variables:
(pn - p0). (3.25)
The resulting system of difference equations is solved by the conjugate gradient method. To accelerate its convergence, the algebraic multigrid method [33] is used.
t
4. Testing of calculation schemes and algorithms
Testing of the correctness of the RANS turbulence models is performed using the problem (having an experimental description [38]) of heat exchange on the inlet surface of the recessed nozzle of the solid propellant rocket engine with a cylindrical channel charge (Fig. 5).
n
Fig. 5. Scheme of the solid propellant rocket engine with a recessed nozzle and a cylindrical charge [38] for testing the RANS turbulence models
According to the experimental data [37, 38], the working substance is air \ M = 0.02896 Cp = 1004.4 p = 1.831 • 10"5 A = 0.0261 the structural material of the model
recessed nozzle is steel (M = 0.05585 p = 7854 ^p-, A = 60.5 The symmetric position
of the recessed nozzle is considered. The influence of the axial displacement of the nozzle unit on gas dynamics is investigated in [5, 6, 37], and the influence of the motion of the recessed nozzle is examined in [45]. The boundary conditions are defined as follows:
On the injection surface (ri), the stagnation pressure, the stagnation temperature and the
sSo'^2' e = 7'
turbulence characteristics are given: p0 = 4.7-105 Pa, T0 = 600 K, k = ^h-U2, s = pC,t ■ —
pt = 1000 • Tu ■ p, Tu = 0.05. On solid impermeable surfaces T3, the adhesion and no-leak conditions ^u = v = w = 0,
du-
= 0 1 are given.
Between the boundary of the solid body and gas (Fig. 2), the boundary condition of the 4th kind = -A, Tw = Tfl) is imposed.
.p
On the nozzle section (r2), soft boundary conditions are given: = 0, f = (p, U).
On the external surfaces of the walls of the combustion chamber (Fig. 2), the boundary condition of the 3rd kind A^f = a(Tw — Ta)j is imposed.
Investigations of mesh independence are conducted in the quasi-stationary mode of operation of the solid propellant rocket engine, i.e., with constant geometry (invariant with respect to the time and flow of gas) of the solid propellant charge channels. Consideration is given to the influence of the computational grid (the ratio of the volume of the largest unit grid element to the general volume of the computational region) on the value of the averaged static pressure in the combustion chamber (Fig. 6).
2-106 4-106 6-106 The number of cells (a)
8-10'
2-10"7 4-10"7 6-10"7 8-10"7 l-lO"6 vjV
(b)
Fig. 6. Graph showing the dependence of the value of the averaged static pressure at the exit from the mass flux channel on (a) the capacity of the grid and (b) the ratio of the volume of the maximal unit grid element to the general volume of the computational region
The curve shown in Fig. 6 demonstrates the stable mesh independence of the problem considered, and the resulting sufficient ratio of the volume of the unit grid element to the general volume of the computational region can be used for discretization with respect to the space of other computational regions.
The computational grid common for all the turbulence models under consideration contains more than 3 million hexahedral cells, including near-wall layers consisting of 82 000 prismatic elements. A computational grid, an enlarged image of the boundary layer in the radial direction, as well as the ratio of the position of the first computational cell relative to the thickness of the boundary layer, which is characterized by the parameter y+, are shown in Fig. 7. All computational schemes, algorithms, and the mathematical model used were verified by considering problems having a detailed experimental description [5, 6].
Figure 8 shows graphs of distribution of the heat-transfer coefficient over the generating line of the recessed nozzle (the coordinate s is laid off from the critical section of the nozzle along its generating line (0 = 00), so that s = 0 is the critical section, s = 1.5 is the top of the nozzle, the region s = 0.5-2 is the inlet of the recessed nozzle), which have been calculated using various turbulence models.
Figure 8 shows that the application of the RNG k—e, k—e and Spalart-Allmaras turbulence models leads, on the one hand, to overestimation of the intensity of heat exchange by 20-50 % in the region of high-speed flow (s ^ 0.5) and to a considerable underestimation (up to 90%) of the values of the heat-transfer coefficient in the region of ingress of the flow to the inlet of the recessed nozzle (0.5 ^ s ^ 2), including the gap above the nozzle (s ^ 1.5) — the zone of reverse flows.
Since the combustion chamber of the solid propellant rocket engine includes connections (gas pipes of different lengths) and traction vector control elements, or short nozzles or the inlet of the recessed nozzle, i.e., regions of high-speed (including trans- and supersonic) flows and
Yplus Contour 1
1
3.830e+000 3.628e+000 3.425e+000 3.223e+000 3.020e+000 2.818e+000 2.615e+000 2.413e+000 2.211e+000 2.008e+000 1.806e+000 1.603e+000 1.401e+000 1.199e+000 9.962e —001 7.937e —001 5.913e —001 3.889e —001 1.865e —001
(b)
Fig. 7. Computational grid (a), prismatic near-wall layers (b) and distribution y+ over the inlet surface of the recessed nozzle (c)
s
Fig. 8. Distribution of the heat-transfer coefficient over the generating line of the recessed nozzle, obtained within the framework of: 1 — the experiment [38]; 2 — the k — w SST model; 3 — the k — w model; 4 — the RNG k — e model; 5 — the k — e model; 6 — the Spalart-Allmaras model; 7 — the SST model
stagnation regions with zones of reverse flows [5], the application of the RNG k — e, k — e and Spalart - Allmaras models for modeling of intrachamber processes is inexpedient.
To check the applicability of the к — ш SST, к — ш and SST turbulence models, an analysis has been made of the problem (having a detailed experimental description [37]) of the flow of combustion products in the flow channels of the solid propellant rocket engine with a recessed hinged nozzle and a starlike solid propellant charge (Fig. 9).
Fig. 9. Scheme of the solid propellant rocket engine with a recessed nozzle and a starlike charge
A computational grid, an enlarged image of the boundary layer, and the ratio of the position of the first computational cell to the thickness of the boundary layer, which is characterized by the parameter y+, are shown in Fig. 10.
(a)
Solver Yplus Contour 1 3.225e+001 3.049e+001 2.873e+001 2.698e+001 2.522e+001 2.346e+001 2.170e+001 1.995e+001 1.819e+001 1.643e+001 1.467e+001 1.291e+001 1.116e+001 9.398e+000 7.640e+000 5.882e+000 4.124e+000 2.366e+000 6.080e —001
(c)
Fig. 10. Computational grid (a), prismatic near-wall layers (b) and the distribution of y+ over the inlet surface of the recessed nozzle (c)
The symmetric position of the recessed nozzle is considered. The ratio of the consumption
of the gas coming from the gap above the nozzle to the consumption of the gas coming from the
G 1
main channel is given by the coefficient kg = jf - = j.
The working substance is air (m = 0.02896 Cp = 1004.4 /x = 1.831 • 10"5
A = 0.0261 ^ j, the structural material of the model recessed nozzle is steel (^M = 0.05585
p = 7854 A = 60.5 The symmetric position of the recessed nozzle is considered. The
influence of the axial displacement of the nozzle unit on gas dynamics is investigated in [5, 6, 37],
and the influence of the rotation of the recessed nozzle is examined in [45]. The boundary conditions are defined as follows:
• On the injection surface (T1), the stagnation pressure, the stagnation temperature and the
turbulence characteristics are given: p0 = 7.8-105 Pa, T0 = 400 K, k = -U2, e = pCj pt = 1000 • Tu ■ p, Tu = 0.05.
On the injection surface in the gap above the nozzle (T2), the gas consumption, the stag-
r0
nation temperature and the turbulence characteristics are given: G = 2-10 Tn = 400 K,
k = m • u'2> 6 = pC» • J' fit = 1000 •Tu • Tu = °-05-
On the solid impermeable surfaces T3 and r5, the adhesion and no-leak conditions (u =
On-
„ du.
= v = w = 0, = 0 J are given.
Between the boundary of the solid body and gas (Fig. 2), the boundary condition of the
4th kind ( = ,
Tw = Tg j is imposed.
On the nozzle section (r5), soft boundary conditions are given:
d2f dn2
0, f = (p, U).
• On the external surfaces of the walls of the combustion chamber (Fig. 2), the boundary condition of the 3rd kind A^f = a(Tw — Ta)j is imposed.
As a result of numerical modeling, spatial distributions of gas-dynamical and thermophysical quantities in the combustion chamber with a starlike charge have been obtained. The structure of the boundary streamlines on the inlet surface of the recessed nozzle, which was obtained by applying various turbulence models, is depicted in Fig. 11.
Fig. 11. The structure of the flow near the inlet surface of the recessed nozzle: (a) experiment [37]; (b) the k - u SST model; (c) the k - u model; (d) the SST model
2
k
Figure 11 shows that, despite the similarity of the structure of the boundary streamlines, the formation of nodes oriented by the rays of the charge compensators, the k — w SST model makes it possible to obtain a topology of gas flow near the inlet surface of the recessed nozzle that is the closest to the experimental soot-oil pictures (one can see the zone of ingress of the flow above the nozzle into the nozzle; this zone is due to a horseshoe vortex surrounding the system of stagnation points).
Figure 12 shows graphs of distribution of the heat-transfer coefficient over the generating line of the recessed nozzle, calculated using various turbulence models. It should be noted that
3
0
2
3
4
5
s
Fig. 12. Distribution of the heat-transfer coefficient over the generating line of the recessed nozzle, obtained within the framework of: 1 — the experiment [37]; 2 — the k — w SST model; 3 — the k — w model; 4 — the SST model
the coordinate s is laid off along the generating line of the recessed nozzle from the point of attachment of the nozzle to the housing of the combustion chamber in the gap above the nozzle
Figure 12 shows that all the turbulence models considered allow one to obtain correct distributions of the heat-transfer coefficient in the region of stagnation points and reverse flows, but, with the model coefficients used, do not correlate with the experimental data in the regions of transsonic and supersonic flows (for s ^ 3.8, Fig. 12 in the critical and expanding parts of the nozzle).
In the subsonic region, the distributions of the heat-transfer coefficient obtained within the framework of the turbulence models under consideration are close to the experimental ones; however, the k — w and SST models give a 16 % error in the determination of the heat-transfer coefficient in the stagnation region near the nozzle bottom (for s ^ 2.2, Fig. 12). Also, the k — w SST model (the 2003 Menter model) in the region of the nozzle bottom provides a satisfactory quantitative agreement with the experiment (with an error of less than 1 %), and the deviation of the calculated data obtained by using the k — w SST model from the experimental data in the region of subsonic flows (for 3.3 ^ s ^ 3.7) does not exceed 8%.
5. Conclusion
An analysis has been made of the applicability of the RANS turbulence models for the modeling of spatial quasi-stationary intrachamber processes in the flow channels and in the pre-nozzle volume of the solid propellant rocket engine. A validation of the most common turbulence models was performed for the problems of conjugate heat exchange using known experimental data. It was shown that the SST k — w turbulence model allows one to obtain local distributions (the closest to the experimental data) of the heat-transfer coefficient in the region of subsonic and reverse flows: the disagreement between the model and the experimental data does not exceed 8 %.
(Fig. 9).
References
[1] Aliev, A. V., Mishchenkova, O. V., and Cherepov, I. V., Nonstationary Intra-Chamber Processes in Solid-Propellant Controlled Propulsion System, Vestn. MGTU. Ser. Mashinostr, 2016, vol. 4(109), pp. 24-39 (Russian).
[2] Bardina, J.E., Huang, P. G., and Coakley, T. J., Turbulence Modeling Validation, Testing, and Development, NASA Technical Memorandum 110446 (Apr 1997).
[3] Bendersky, B. Ya. and Teneev, V. A., Spatial Subsonic Flows in Areas with Composite Geometry, Matem. Mod., 2001, vol. 13, no. 8, pp. 121-127 (Russian).
[4] Bendersky, B.Ya., Cherepov, I.V., and Il'yashenko, K.V., The Study of Gas Dynamic Processes in Engine Management, Vestn. MGTU. Ser. Mashinostr., 2004, no. 5, pp. 95-103 (Russian).
[5] Benderskiy, B.Ya. and Chernova, A.A., Formation of Vortex Structures in Channels with Mass Injection and Their Interaction with Surfaces in Solid-Fuel Rocket Engines, Thermophys. Aeromechanics, 2015, vol. 22, no. 2, pp. 185-190.
[6] Benderskiy, B. Ya. and Chernova, A. A., Features of Heat Transfer in a Pre-Nozzle Volume of a Solid-Propellant Rocket Motor with Charges of Complex Shapes, Thermophys. Aeromechanics, 2018, vol. 25, no. 2, pp. 265-272.
[7] Bendersky, B. Ya. and Chernova, A. A., Investigation of Heat Transfer in the Combustion Chamber of a Solid-Propellant Rocket Motor within the Framework of the Homogeneous Gas Model, Khimich. Fiz. Mezoskop, 2021, vol. 23, no. 4, pp. 412-429 (Russian).
[8] Garbaruk, A. V., Current Approaches to Turbulence Modeling, St. Petersburg: SPbPU, 2016 (Russian).
[9] Garbaruk, A. V., Viscous Fluid Flows and Turbulence Models: Methods for Calculating Turbulent Flows, St. Petersburg: SPbPU, 2010 (Russian).
[10] Garbaruk, A., Strelets, M. Kh., and Shur, M. L., Turbulence Modeling in Complex Flow Calculations, St. Petersburg: SPbPU, 2012 (Russian).
[11] Isaev, S., Gritckevich, M., Leontiev, A., and Popov, I., Abnormal Enhancement of Separated Turbulent Air Flow and Heat Transfer in Inclined Single-Row Oval-Trench Dimples at the Narrow Channel Wall, Acta Astronaut, 2019, vol. 163, part A, pp. 202-207.
[12] Isaev, A. I. and Skorobogatov, S. V., Hydrodynamic Verification and Validation of Numerical Methods of the Flow Calculation in Combustion Chamber of a Gas Turbine Engine, Trudy MAI, 2017, no. 97, 28 pp. (Russian).
[13] Jagadeesh, P. and Murali, K., Application of Low-Re Turbulence Models for Flow Simulations Past Underwater Vehicle Hull Forms, Int. J. Nav. Archit. Ocean Eng., 2005, vol. 2, no. 1, pp. 41-54.
[14] Kolmogorov, A. N., Equations of Turbulent Motion of an Incompressible Fluid, Dokl. Akad. Nauk SSSR. Ser. Fiz., 1942, vol. 6, nos. 1-2, pp. 56-58 (Russian).
[15] Koroleva, M.R., Mishchenkova, O.V., Kelemen, M., and Chernova, A. A., Theoretical Research of the Internal Gas Dynamics Processes of Measurements of Hot Air Curtain with Cross-Flow Fan, MM Sci. J., 2020, June, pp. 3966-3972.
[16] Koroleva, M.R., Mishenkova, O.V., Raeder, T., Tenenev, V. A., and Chernova, A. A., Numerical Simulation of the Process of Activation of the Safety Valve, Kompyuternye Issledovaniya i Mod-elirovanie, 2018, vol. 10, no. 4, pp. 495-509.
[17] Gas-Dynamic and Thermophysical Processes in Solid Propellant Rocket Engines, A. S.Koroteev (Ed.), Moscow: Mashinostroenie, 2004 (Russian).
[18] Kravchuk, M.O., Kudimov, N.F., and Safronov, A.V., Turbulence Modeling Issues for Supersonic High-Temperature Jet Computation, Trudy MAI, 2015, no. 82, 20 pp. (Russian).
[19] Larina, E. V., Kryukov, I. A., and Ivanov, I. E., Numerical Simulation of Axisymmetric Jets Using Differential Eddy Viscosity Models, Trudy MAI, 2016, no. 91, 24 pp. (Russian).
[20] Launder, B.E. and Spalding, D.B., Lectures in Mathematical Models of Turbulence, London: Acad. Press, 1972.
[21] Internal Ballistics of Solid Propellant Rocket Engines, A. M.Lipanov et al., Moskow: Mashinostroe-nie, 2007 (Russian).
[22] Lipanov, A.M., Bobryshev, V. P., Aliev, A. V., Spiridonov, F.F., and Lisitsa, V.D., Numerical Experiment in the Theory of Solid Propellant Rocet Engines, Ekaterinburg: Nauka, 1994 (Russian).
[23] Lipanov, A.M., Dadikina, S.Yu., Shumikhin, A.A., Koroleva, M.R., and Karpov, A.I., Numerical Simulation Intra-Chamber of Unsteady Turbulent Flows Stimulate: Part 1, Vestn. YuUrGU. Ser. Mat. Model. Progr, 2019, vol. 12, no. 1, pp. 32-43 (Russian).
[24] Menter, F.R., Zonal Two Equation k — w Turbulence Models for Aerodynamic Flows, in Proc. of the 24th Fluid Dynamics Conf. (American Institute of Aeronautics and Astronautics, Orlando, FL, July 1993), pp. 128-143.
[25] Menter, F.R., Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications, AIAA J, 1994, vol. 32, no. 8, pp. 1598-1605.
[26] Menter, F.R. and Esch, T., Advanced Turbulence Modelling in CFX, CFX Update, 2001, no. 20, pp. 4-5.
[27] Menter, F. R., Kuntz M., and Langtry, R., Ten Years of Industrial Experience with the SST Turbulence Model, in Turbulence, Heat and Mass Transfer: Proc. of the 4th Internat. Symp. on Turbulence, Heat and Mass Transfer (Antalya, Turkey, Oct 2003), K. Hanjalic, Y. Nagano, M. J. Tummers (Eds.), Danbury, Conn.: Begell House, 2003, pp. 625-632.
[28] Grotjans, H. and Menter, F., Wall Functions for General Application CFD Codes, in ECCOMAS'98: Proc. of the 4th European Computational Fluid Dynamics Conf., K. D. Papailiou et al. (Eds.), Chichester: Wiley, 1998, pp. 1112-1117.
[29] Naplekov, I. S., Computer Modeling for Solving the Problem of Gas Flow from a Nozzle Using Different Turbulence Models, in The 14th Korolev's Readings: The Internat. Young Scientists Conf. (Samara, Oct 2017), Samara: SGU, 2017, pp. 209-210 (Russian).
[30] Orszag, S.A., Yakhot, V., Flannery, W. S., Boysan, F., Choudhury, D., Maruzewski, J., and Pa-tel, B., Renormalization Group Modeling and Turbulence Simulations, in Proc. of the Internat. Conf. on Near- Wall Turbulent Flows (Tempe, AZ, March 1993), R. M. C. So, Ch. G. Speziale, B. E. Launder (Eds.), Amsterdam: Elsevier, 1993, pp. 1031-1046.
[31] Raeder, T., Tenenev, V. A., Koroleva, M. R., and Mishchenkova, O. V., Nonlinear Processes in Safety Systems for Substances with Parameters Close to a Critical State, Russian J. Nonlinear Dyn., 2021, vol. 17, no. 1, pp. 119-138.
[32] Raeder, T., Tenenev, V. A., and Chernova, A. A., Incorporation of Fluid Compressibility into the Calculation of the Stationary Mode of Operation of a Hydraulic Device at High Fluid Pressures, Russian J. Nonlinear Dyn., 2021, vol. 17, no. 2, pp. 195-209.
[33] Trottenberg, U., Oosterlee, C. W., and Schüller, A., Multigrid, New York: Acad. Press, 2000.
[34] Ferziger, J.H. and Peric, M., Computational Methods for Fluid Dynamics, 2nd rev. ed., Berlin: Springer, 1999.
[35] Safronov, A.V., Applicability of Turbulent Viscosity Model for the Calculation of Supersonic Jet Streams, Fiz.-Khim. Kinetika Gaz. Dinam., 2012, vol. 13, no. 1, Art. 9, 17 pp. (Russian).
[36] Sandhem, N. D. and Kleiser, L., The Late Stages of Transition to Turbulence in Channel Flow, J. Fluid Mech., 1992, vol. 245, pp. 319-348.
[37] Saveliev, S.K., Emelyanov, V. N., and Benderskiy, B.Ya., Experimental Methods of SRM Gas Dynamics Research, St. Petersburg: Nedra, 2007 (Russian).
[38] Shishkov, A. A., Panin, S. D., and Rumyantsev, V. V., Workflows in Solid Propellant Rocket Engines, Moskow: Mashinostroenie, 1989 (Russian).
[39] Spalart, P. R. and Allmaras, S.R., A One-Equation Turbulence Model for Aerodynamic Flows, Technical Report AIAA-92-0439, American Institute of Aeronautics and Astronautics, 1992, 22 pp.
[40] Starchenko, A. V., Nuterman, R. B., and Danilkin, E. A., Numerical Modeling of Turbulent Currents and Impurity Transport in Street Canyons, Tomsk: TGU, 2015 (Russian).
[41] Launder, B.E. and Spalding, D.B., The Numerical Computation of Turbulent Flows, Comput. Methods Appl. Mech. Eng., 1974, vol. 3, no. 2, pp. 269-289.
[42] Volkov, K. N. and Emelyanov, V. N., Mass-Driven Gas Flows in Channels and Ducts of Power Plants, Moscow: Fizmatlit, 2011 (Russian).
[43] Volkov, K. N. and Emelyanov, V. N., Mathematical Models of Three-Dimensional Turbulent Flows in the Ducts with Fluid Injection, Matem. Mod., 2004, vol. 16, no. 10, pp. 41-63 (Russian).
[44] Volkov, K. N. and Emelyanov, V. N., Large Eddies Simulation in Turbulent Flow Calculations, Moscow: Fizmatlit, 2011 (Russian).
[45] Volkov, K. N., Denisikhin, Emelyanov, V. N., and Teterina, I. V., Numerical Simulation of Gas-Dynamics Processes in Thrust Vectorable Nozzle, Fiz.-Khim. Kinetika Gaz. Dinam., 2018, vol. 19, no. 2, 24 pp.
[46] Wilcox, D.C., Turbulence Modeling for CFD, La Canada, Calif.: DCW Industries, 1998.
