CC BY
4Russian Journal of Nonlinear Dynamics, 2021, vol. 17, no. 1, pp. 119-138. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd210109
NONLINEAR ENGINEERING AND ROBOTICS
MSC 2010: 35Q30, 76-10, 76N15
Nonlinear Processes in Safety Systems for Substances with Parameters Close to a Critical State
T. Raeder, V. A. Tenenev, M. R. Koroleva, O. V. Mishchenkova
The paper presents a modification of the digital method by S.K.Godunov for calculating real gas flows under conditions close to a critical state. The method is generalized to the case of the Van der Waals equation of state using the local approximation algorithm. Test calculations of flows in a shock tube have shown the validity of this approach for the mathematical description of gas-dynamic processes in real gases with shock waves and contact discontinuity both in areas with classical and nonclassical behavior patterns. The modified digital scheme by Godunov with local approximation of the Van der Waals equation by a two-term equation of state was used for simulating a spatial flow of real gas based on Navier-Stokes equations in the area of a complex shape, which is characteristic of the internal space of a safety valve. We have demonstrated that, under near-critical conditions, areas of nonclassical gas behavior may appear, which affects the nature of flows. We have studied nonlinear processes in a safety valve arising from the movement of the shut-off element, which are also determined by the device design features and the gas flow conditions.
Keywords: real gas, Van der Waals equation, critical state of substance, Godunov's method
Received March 02, 2021 Accepted March 15, 2021
The work was performed as part the research project of the Udmurt Federal Research Center of the Ural Branch of the Russian Academy of Sciences "Artificial Intelligence in the Development, Training, and Support of Expert Systems for Presentation and Use of Knowledge in Natural, Engineering, and Socio-Humanitarian Science" AAAA-A19-119092690104-4.
Thomas Raeder [email protected] Olga V. Mishchenkova [email protected]
Kalashnikov Izhevsk State Technical University ul. Studencheskaya 7, Izhevsk, 426069 Russia
Valentin A. Tenenev [email protected]
Kalashnikov Izhevsk State Technical University ul. Studencheskaya 7, Izhevsk, 426069 Russia Udmurt Federal Research Center UB RAS ul. T. Baramzinoi 34, Izhevsk, 426067 Russia
1. Introduction
Production processes are directly associated with the use of various liquid and gaseous substances, the storage and transportation of which are often accompanied by emergency situations. This leads to equipment failure or damage. For this purpose, the use of safety systems is a necessary measure to prevent emergency operational conditions. The considered type of safety devices functions as a result of contact with a working medium. The latter can be gas, steam, water, lubricating and hydraulic oils, petroleum products, etc. The type of working medium determines the system features and regulation parameters. Therefore, investigating the medium behavior in the device volume is of great practical importance, especially for the range of changes in the working medium characteristics that covers the critical state.
The Mendeleev-Clapeyron equation of state for ideal gas is widely used in the mathematical description of gas-dynamic processes in safety devices. This equation well describes gas within a wide range of temperature and pressure changes when the medium is under conditions far from condensation. However, when operating conditions do not coincide with the applicability conditions for the ideal gas model, we need to choose models that are correct under conditions close to the critical state of substances. Thus, a method for calculating the thermophysical properties of substances based on the interpolation of tabular data is presented in [29]. The digital simulation of compression and energy release processes in a target with combined effects of a pulsed jet system and intense laser radiation using the method of magnetic-inertial plasma confinement involves calculated dependences of thermodynamic characteristics that determine the equation of state of the substance in a form different from that for ideal gas [15]. Nevertheless, the most commonly used are equations of state for real gas, such as the Van der Waals, Redlich-Kwong, Peng-Robinson, Kamerlingh Onnes, Martin Howe equations, etc. [2, 3, 17, 20, 23, 30]. The best-known equation of state for real gas is the Van der Waals equation, which takes into account intermolecular interactions and the real (nonzero) intrinsic volume of molecules. This model allows for reflecting a real set of different states of gas/liquid, depending on the distance from the critical point. Thus, based on the Van der Waals polytropic gas model, nonclassical gas-dynamic behavior of vapors near the liquid-vapor saturation curve is studied in [5] and a modification of the Van der Waals equation is used to describe explosives within the range of pressures and temperatures characteristic of detonation processes [13]. Based on the digital solution to a system of Navier-Stokes equations with the closure of the system by the Van der Waals equation of state, theoretical studies of flows and heat transfer in near-critical media under space flight conditions are carried out in [21]. A modified Van der Waals gas model is used to analyze the mass flow rate and opening parameters of safety valves at operating pressures up to 3,600 bar [2], and to study flows in engineering devices used to implement the organic Rankine cycle [20].
In this paper, we present a digital simulation of the gas flow in a direct-acting safety valve using the Van der Waals equation of state. The working medium is C13F22 solvent used in selective CO2 capture processes from fuel gas flows [7]. The specifics of gas-dynamic process digital simulation are in the thermophysical properties of high-molecular substances. Near the critical point of the substance, its isentropes lose their convexity; as a result, there is an increased probability of appearance of nonclassical flow regions.
Maria R. Koroleva [email protected]
Udmurt Federal Research Center UB RAS ul. T. Baramzinoi 34, Izhevsk, 426067 Russia
Numerical simulation is based on the S. K. Godunov's method. Generalization of the Godunov-type schemes [9] to the case of general equation of state faces several problems. The main one is the rejection of numerous references to the equations of state in the process of solving the Riemann problem and minimizing the used information about the equation of state. For example, an algorithm similar to that described in [19], which requires isentrope integration, can lead to significant computational costs, on the one hand, and to the loss of the method stability in depression regions, on the other hand. To overcome these difficulties, in [31] the sound speed in the equation of state was obtained by digitally differentiating the tabular pressure values, which provided the same stability and reliability in calculations. A modification of Rowe's linearization method for real gases and its application to a particular case of the Van der Waals gas are presented in [10] and presume introducing an intermediate state as an unknown value of the linearization problem. Due to the analytical form of the Van der Waals equation of state, we have obtained a third-order algebraic equation for the intermediate density, which directly provides a solution to the linearization problem in terms of Rowe mean velocity and total enthalpy, as well as jumps in density, pressure, and internal energy per unit volume. The paper [8] considers the system of Euler equations for the Van der Waals polytropic gas. The solution to the Riemann problem is reduced to a system of two nonlinear equations for finding the density of real gas on both sides of the contact discontinuity, giving solutions to the Riemann problem for three cases, including nonclassical and mixed waves.
In this paper, we present a modification of the digital method by S.K.Godunov for calculating the flows of real gases under conditions close to a critical state. The method is generalized to the case of the Van der Waals equation of state using the local approximation algorithm.
2. Van der Waals equation of state
The Van der Waals equation of state, which describes the relationship between the specific energy e, entropy S and density p of a substance, is as follows [23]:
e (5, p) = -Ap + (eo + Apo) (\~_Bßp ) (£) exp((S-So)d/R), (2.1)
where A and B are dimensional Van der Waals constants, depending on a specific substance, d = R/Cv, R is the specific gas constant, Cv is the specific heat at constant volume; constants with a 0 index correspond to some initial (original) state of the gas.
From Eq. (2.1) we can obtain an expression for the absolute pressure
P (S, p) = (?de^p) = -Ap2 + d (so + Ap0) {1~Bpo)d pl+d ((S - S0) d/R).
dp pd (1 - Bp)1+d
After excluding the entropy from the obtained formula, using relation (2.1), we obtain a dependence of the pressure on the gas density and specific energy
p(e,p) = pd^±^-Ap2. (2.2)
1 - Bp
This equation can be represented in terms of the medium's temperature T as follows:
d
T(s,p) = -(s + Ap). (2.3)
In the case of an isentropic process, Eqs. (2.1) and (2.2) can be represented as
e (p) = -Ap + (eo + Apo)
' 1 - Bpo
l-BppY 1 -Bp)
P (p) = (po + Api)
1 - Bp
1+d
P_ p0
1+d
- Ap2
When calculating gas flows, it is convenient to represent the initial equations in dimensionless form by means of the ratio of gas-dynamic quantities to some characteristic values. In the case of real gas, it is preferable to use the Van der Waals equation of state in dimensionless form with respect to the critical parameters of the substance: pc, pc, Tc.
The dimensional coefficients of the Van der Waals equation
_ 27 R2T2
~ 64 pc '
8 Pc
after the nondimensionalization procedure, will be equal to: a = 3, b = 1/3 [5], i.e., they do not depend on the medium in question. Then the equation for the internal energy, which depends on temperature and density, takes the dimensionless form
e (T,p) =
8bT
ap.
(2.4)
3. Mathematical simulation of a real gas
in an arbitrary discontinuity disintegration problem
The breakdown of an arbitrary gas-dynamic discontinuity problem [9], or the Riemann problem, is related to solving a system of Euler one-dimensional gas-dynamic equations:
dQ_ dF_ dt dx
0,
(3.1)
p pu
Q = pu , F = pu2 + p
e (e + p) u
where t is time, x is a coordinate, the vectors Q and F contain gas-dynamic complexes composed of variables, p is the gas density, u is velocity, e = p(e + u2/2) is the total energy per gas unit volume, and e = e(p, p) is the internal energy per unit mass given by the general equation of state.
The values of gas-dynamic complexes at the initial moment of time are set in the regions to the left and right of the discontinuity:
Q (0, x) =
Qi, x < 0, Qi, x> 0.
(3.2)
The flow structure is determined by a combination of either the emerging rarefaction waves and shock waves with a contact discontinuity, or the vacuum region.
d
d
The flow characteristics are fully described by the well-known relations for the shock wave and the rarefaction wave. For a shock wave, the jump of parameters is determined by the Rankine-Hugoniot conditions [4]:
l»] - ±M,
m
[p]
m2
M = -72(P+P)
(3.3)
(3.4)
(3.5)
where [u] = U — u, [p] = P — p, [e] = E — e, [1/p] = 1/R — 1/p. Here, u, p, e, p are velocity, pressure, internal energy, and density, respectively, ahead of the shock wave front, whereas U, P, E, R are the same values behind the shock wave front, and m is the mass velocity. The following equations are valid for the depression wave:
dp dp = c2, (3.6)
du dp 1 = ±—, pc (3.7)
de dp P {p, e) P2 ' (3.8)
where c is the sound speed.
In addition to relations (3.6)-(3.8), we use the continuity condition for the contact discontinuity of velocity U and pressure P. From conditions (3.3)-(3.5), we have obtained a system of four equations for determining four unknown values on the contact discontinuity U, P, R1, R2. Relations (3.3) give us
U — u1 = —
P — Pi
U — U2 =
P — P2
mi m2
where we obtain expressions for the mass velocity from Eqs. (3.4) and (3.7)
(3.9)
'RiPi^-^-, P>pu
Ri — pi'
mi =
P-Vi P
r CJ£
J pc
P < Pi,
(3.10)
Pi
corresponding to (3.3) for the shock wave P > p or the Poisson adiabat (3.8) for rarefaction waves P < p. Index i = 1,2 in Eq. (3.10) corresponds to one of the initial states of gas according to the initial conditions (3.2).
After excluding the velocity U from relations (3.9), we obtain a nonlinear equation connecting the parameters P, R1, R2 :
/i (P, Ei, R2) = P ( — + — 1 m1 m2
Pi m1
P 2 m2
— u1 + u2 = 0.
(3.11)
For shock waves, the energy equation is as follows:
E{P,R) + -2{P + p) (1.
where the function E(P, R) corresponds to the equation of state of the substance under consideration.
Generally, in the case of depression waves, we shall digitally integrate the differential equation of the Poisson adiabat (3.8).
The second and third equations for unknown values of P, Ri, R2 are as follows:
/2 (P, Ri) = 0 =
/3 (P, R2) = 0 =
P - Ps (pi,pi,Ri),
e(p,r.2) + \(p + p2) p2 p - Ps (P2,P2,R2),
1
Pi 1
P>Pi, P <Pi, P>P2, P <P2,
where ps is the pressure in the rarefaction wave.
As a result, we obtain a system of nonlinear equations:
/i(P,Ri,R) = 0, /2(P, Ri) = 0,
h(P, R2) = 0,
(3.12)
which is solved at some initial approximation of (P, R1, R2)0 using the Newton method. Calculating the mass velocity m according to (3.10) in the case of rarefaction wave requires calculating
p
the integral j —, which can be represented in quadratures only for some special cases [18]. p
The value of speed U after calculating the values of P, R1, R2 can be found from the conditions
rr P-Pl . P-P'2
U = II, 1--= 112 +
mi
m2
The exact solution to the Riemann problem in the case of equations of state for real gas is associated with an iterative solution to the system of equations (3.12), which is an expensive operation, in terms of computational costs, in the digital implementation of Godunov-type difference schemes. There is a known approach with local approximation of an equation of state by a two-term equation [14, 18, 22], where relation (2.2) can be reduced to [9]:
£ (p, p) =
p - 4 (p - pk) (7-1 )P '
(3.13)
where 7 is the adiabatic exponent and ck , pk are constant values. Equation (3.13) satisfies the restrictive Bethe-Weil inequalities [9] and can be used in methods based on solving the problem of arbitrary initial discontinuity disintegration, which include the classical Godunov method. For a two-term equation of state, the exact Riemann solver was described in detail by S. K. Godunov,
being used as the basis for the difference schemes. In this case, the Van der Waals equation (2.2) is represented as follows [28]:
e(p, p) = + Ap^j - Ap. (3.14)
The corresponding expressions for the change in pressure and sound speed in the rarefaction wave take the form
- for the two-term equation of state (3.13) [9]:
a (S) RY pkc\ p + pk
ps{p,p,R) =--pk, pk =-<7(5) =7——,
Y Y pY
, p + pk
for the Van der Waals equation (2.2) [23]:
p. (p, p, R) = {p + Ap2) (I)" M2
c = ^(l + d)[j-2+A)T^--2Ap.
For the two-term equation of state, we obtained an analytical expression for the mass velocity for the rarefaction wave [9]:
1 _ P+Pfc 7-1 P + Pk
m = -^TPc-Hi'
1 _ (P + Pk\ 27
\P +Pk
which eliminates the need for costly digital integration of equations (3.11).
If the equation of state is given as p(e, p) = 0, it can be represented as p (e, p) =
= t — —pc% parame^-ers ^ Cfc iocany determined from conditions [14]:
P (Y — 1)
= i =--P2( 7-1) =°' ^ =
Local approximation of parameters y, po, co is
, . 1 dp 1 i dp \ 2 p (e,P)+ pk f -.x
7 = 1 + -^, Pk = ~ [ Pt^-P{£,P) ), ck =---7T--t 7~l -
pdp Y \ dp J p (y — 1)
Reducing the Van der Waals equation to a two-term equation of state, we can obtain the following expressions:
d 1 ((1 + d) Ap2 + p „ 2 7 = 1 + ^-J5-, Pk = - --t2—^-- - 2 Ap2 - p
1 - Bp 1 7 V 1 - Bp
A local approximation of the equation of state allows us to find the initial approximation of variables (P, R1, R2)0 to solve the system of equations (3.12). When performing test calculations, the system of equations was solved until the disparity |A| ^ 10-7 was reached. To include the algorithm for solving the Riemann problem for real gases in the difference method due to S.K.Godunov, we need to check the sufficiency of the initial approximation. The results of numerous test calculations have confirmed the possibility of using the initial approximation based on the local parameter approximation of a two-term equation of state, which can also be used in a tabular setting of equations of state for substances. The application scope of the local approximation is determined by the following conditions:
dp „dp d2p -f>0, 2 -rf+p-dp dp t
where T is the gas temperature.
dp dS
-^■>0, 2^+p—1>0, 7^>0, — >0, lim p (S, p) —>■ oo, dp dp op2 -
de dT
p^-IX
4. Test calculations for real gas
All calculations have been performed in dimensionless form. The nondimensionalization procedure has been carried out based on the critical parameters of the substance: pc, pc, Tc. In this case, internal energy should be attributed to the value pc/pc, and the gas velocity, to sj'pdpc- The Van der Waals equation of state in dimensionless form is represented by Eq. (2.4).
Test calculations have been performed for the problem of flow in a shock tube with initial parameters in the left half with x < 0: p1, p1, u1 and in the right half with x > 0: p2, p2, u2 for the Van der Waals equation of state. The range of changes in the flow characteristics includes the near-critical area of substance with d = 0.329. Figure 1 shows a comparison of the digital calculation results with the exact solution to the arbitrary discontinuity disintegration problem [9] in the classical field. The lines show the exact solution, the markers, the digital one. At the initial moment of time, the gas parameters were set as follows: p1 = 1, u1 = 1, p1 = 0.4 and p2 = 1, u2 = —0.5, p2 = 0.8. Following the discontinuity disintegration, two oppositely directed shock waves are formed in gas, and the contact discontinuity boundary is shifted to the right relative to the initial position. A comparison of the calculated parameters U, P, R1, R2 with the results of the study [23] is shown in Table 1. One can see the correspondence of the calculated parameters from the solution of Eqs. (3.12) to the values in [23] and the exact solution in [9]. The maximum difference between the results obtained using the local approximation method and the exact solution was 0.002.
To study the nonclassical behavior of the Van der Waals gas near the critical point, we have considered test problems from [1, 10]. The nonclassical behavior of gas is characterized
Table 1. Comparison of parameters U, P, R1, R2
P U Ri Ra
exact solution 1.91913 0.0775 0.63527 1.12721
local approximation 1.916914 0.07606 0.634292 1.126308
results [23] 1.91913 0.077503 0.63527 1.1272
2 1.5
S
PM
0.5 0
-0.15
(a)
0
x
"1
0.15
2 1.5
» 1 s 1
u
Q
0.5 0
(b)
-0.15
0
x
0.15
1.5 1
£ 0.5
u
0
1 0
-0.5 -1
-0.15
(c)
0
x
0.15
Fig. 1. Comparison of exact and digital solutions to the problem with two shock waves, pressure (a), density (b), velocity (c); the solid line marks the exact solution, the dots denote the numerical solution.
by a negative value of the fundamental derivative [12]:
G(v,p) = -,v
~d2p(S,v)
du2 S
[dp(S»l
du S
where v is the specific volume which characterizes the isentrope curvature in the plane (v,p).
' d2p (S,v)'
If the strict inequality
du2
> 0 is satisfied, the isentrope is a convex function. The
isentrope is a straight line if G = 0, and a concave curve if G < 0.
Classical gas dynamics is based on the assumption that the fundamental derivative is positive, and an example thereof is the formation of a shock wave by merging plane compression waves of a finite amplitude. On the contrary, nonclassical gas-dynamic phenomena are expected in flows developing in the negative (G < 0) or mixed nonlinearity, where the sign of the derivative changes from negative to positive and vice versa. Nonclassical phenomena that arise in this case in the medium include the formation of expansion waves, shock waves upstream or downstream of the flow, as well as mixed and separated waves [11].
For the Van der Waals equation, this characteristic is presented as [10]:
G (p,p) =
2
(1 + d) (2 + d) p2 p + apr - 6apA
_(1 -bp)2
2
In our test calculations, we have considered a substance with a high specific heat capacity, relative to its molecular weight, with coefficient d = 0.0128. The calculations were performed for a mixed problem. Figure 2 shows the results of calculating the flow in a shock tube: graphs of changes in the gas pressure, velocity, and density, as well as the fundamental derivative value. The initial state of the gas in this example is characterized by a negative left value Gi = —0.018 and a positive right value G2 = 0.704. The sign change occurs when x ~ —0.1. The gas parameters in the shock tube at the initial moment of time are as follows: p1 = 1.09, u1 =0, p1 = 0.879
x
Fig. 2. Numerical solution to the Riemann problem for a mixed problem, 1 — pressure, 2 — density, 3 — velocity, 4 — fundamental derivative.
and p2 = 0.575, u2 = 0, p2 = 0.275. There is a rarefaction wave arising in gas. However, unlike the classical wave which is formed at positive values of the fundamental derivative, in this case, a mixed wave appears, which is reflected in the profile steepness of all parameters [6]. Part of the front corresponding to negative values of G (x < —0.111) is characterized by a greater steepness, which significantly decreases upon passing to the classical flow region (x > —0.111).
Our test calculations have shown that the use of local approximation in the Godunov method is a good replacement for the exact Riemann solver for real gases both in the classical and nonclassical regions of gas behavior. This approach has been generalized to the three-dimensional case of a real gas flow in an area of complex shape, which is characteristic of the internal space of a safety valve.
5. Gas dynamics of real gas in a safety valve using the Godunov scheme
The diagram of a direct-acting safety valve is shown in Figure 3a. The operation principle of the safety system is to release excess pressure into the environment through the outlet pipe when a certain critical value of the pressure in the inlet channel is reached. The dynamic process occurring in the system under consideration is determined by the thermodynamic parameters of the substance in the reservoir, the characteristics of the spring material, and the spring pretension. The disk movement is determined by both gas dynamics and spring response. A gradual drop in pressure in the reservoir causes the valve to close under the effect of the spring's elastic force [24, 26].
Fig. 3. Direct-acting safety valve: (a) general view: 1 — inlet channel, 2 — outlet pipe, 3 — disk, 4 — spring, (b) calculated region.
The behavior of real gas in a safety valve is studied in a three-dimensional setting. The valve design allows for using an axisymmetric coordinate system for the internal volume of the body. The outlet pipe is calculated in a Cartesian coordinate system.
The mathematical model in the form of a system of differential equations describing the variation of medium parameters in the valve internal volume in dimensionless form for the Cartesian coordinate system is written as follows:
Q =
dQ , d(F - Fv ) _ d(G - Gv ) , d(H - Hv)
Fv =
dt
+
dx
+
dy
+
dz
= 0,
p pVx pVy pVz
pVx p + pvx pVy Vx pVz Vx
pVy , F = pvx vy , G = 2 p + V , H = pVz Vy
pVz pVx Vz pVy Vz P + pV2
e (e + p) Vx _ (e + P) Vy _ (e + p) Vz
0 0 0
Txx Tyx Tzx
Txy 5 Gv = Tyy 5 Hv = Tzy
Txz Tyz Tzz
T xV qx _ TyV - Qy _ rzv - qz
(5.1)
where x, y, z are coordinates, t is time, v = (vx,vy,vz) is the velocity vector, r is the stress tensor, r is the stress vectors, and q is the heat flow.
In a cylindrical coordinate system, the equations are written as follows:
drQ dr(F - Fv) dr(G - Gv) dr(H - Hv)
dt
+
Q =
P
pVx pVr pVg
e
F =
dx
PVx
p + pvX
PVxVr
pVxVe
(e + p) v,
+
dr
G =
+
PVr PVrVx
p + pv;
pVr ve (e + p) v
d9
L,
(5.2)
H
pve
pvd vx pve vr
2
p + pv; (e + p) ve
Fv =
0 0 0 0
Txx Txr Txe 0
Trx , Gv = Trr , Hv = Tre , L = p + pv; - Tee
T0x T$r Tee -pvr ve + Tre
T xv qx trv — qr tqv - qe 0
where r, 0 and vr, ve are radial and angular coordinates and velocity components, respectively. These gas dynamic equations are supplemented by the equations of the k-e turbulence model. A more detailed mathematical setting can be found in [16, 25].
On solid surfaces, we use the no-slip condition; on the exit boundary, we extrapolate the parameters. At the input boundary, we set the stagnation enthalpy and the constant entropy condition. To determine the gas-dynamic parameters at the inlet boundary, we solve the following equations:
- the enthalpy constancy equation:
< \ I p I v2 , x , poo
t (P, P) + - + TT = t {poo, Poo) +-;
P 2 Poo
isentropic pressure equation:
p = ps (poo, Poo, p) = (poo + ap0o)
the Riemann invariants constancy equation:
P
1 ~ bpoo p
1 - bp poo
1+d
ap
vn = ci--(ci - (vn) 1),
Pi
where vn is the speed normal to the boundary, poo and poo are stagnation parameters, and index 1 corresponds to the first boundary cell. Three equations are combined into one: ^(p) = 0, which is solved by the Newton method.
At the initial moment of time, in part of the calculated region in front of the narrow gap, we set the maximum gas pressure and density p = 1 and p = 1, and behind the gap, we set the pressure and density corresponding to environmental parameters p = 0.2 and p = 0.1, and gas is assumed to be stationary. We consider a high-molecular substance — perhydrofluorene Ci3F22 [5, 7] with the following characteristics: pc = 15.8 bar, pc = 1380 kg/m3, Tc = 632 K, d = 0.0128. The calculation has been performed for the minimum and maximum disk lifting heights: 1 mm and 14 mm.
For a numerical solution to the systems of equations (5.1) and (5.2), we have used the control volume method. For each face, we have solved the corresponding discontinuity disintegration problem with two sets of gas-dynamic parameters from the control volumes adjacent to such face. The cylindrical part of the valve and the outlet pipe are joined on the outer generatrix of the cylindrical area. The flows passing through the common faces of control volumes have been calculated using the velocity component normal to the face and thermodynamic values from adjacent control volumes.
The calculation mesh is shown in Fig. 3b with 4-fold rarefaction, since coordinate lines merge on a fine grid. In the cylindrical part, the grid size is 688 x 128 x 48 meshes; in the outlet, the grid is unstructured and contains 354,816 meshes. To increase the order of approximation of the Godunov difference method, we used the MUSCL scheme. The digital solution method is described in detail in [25].
The calculation results are shown in Figs. 4-12.
Exceeding the threshold value of the pressure in the inlet channel rises the disk and releases gas into the valve interior space. The medium is discharged in only one direction: through the outlet pipe. The higher the pressure in the inlet, the higher the force acting on the disk and, accordingly, the greater the valve opening. Figure 4b shows the gas flow structure in the form of stream tubes when the disk rises to a height of 1 and 14 mm. The flow structures in both cases are similar due to the valve's internal cavity geometry. The gas flow entering the valve's closed area splits into two flows after colliding with the body's back wall. You can see that a large-scale recirculation flow area is formed in the lower part. A similar structure of a more complex shape is observed in the upper part above the disk. Earlier, in [27] it was noted that such a flow structure has an additional force effect on the disk and can affect the operating mode of the valve in general. Vortex currents can be observed near the top cover of the body and directly under the disk at the separation of flows. The flow in the exhaust channel contains no large-scale vortices.
-0.1
1 1.2
(a)
(b)
Fig. 4. Flow structure in the valve when opening by: (a) 1 mm, (b) 14 mm.
The color shade of flow tubes corresponds to the fundamental derivative value. If the valve is opened by 1 mm (Fig. 4a), negative values of G(p,p) are observed throughout the inlet channel, where the gas is under high pressure. After the critical section upon passing the sound speed, the derivative becomes positive: a mixed nonlinearity mode is established. The mixed nonlinearity mode is also typical for valve flow with a maximum opening of 14 mm (Fig. 4b). In this case, the derivative changes its sign upstream of the flow, at the level of geometric narrowing of the input channel.
A more detailed representation of changes in gas parameters along the flow is given in Figs. 5-12. Figure 5 shows the change in pressure and flow rate at a minimum — 1 mm (dashed line) — and maximum — 14 mm (solid line) — valve opening. In the figures, the value L is plotted along the lower axis, marking the distance from the valve inlet section along the longitudinal coordinate grid line, corresponding to the boundary layer limit. The critical section is at the point L = 164 mm. One can see that the pressure release occurs faster at a disk lift height of 14 mm: the disturbances propagate much further downstream and upstream of the
(a)
(b)
1
\
2 i i
i \
\ i. / ^ / f
100 200 L (mm)
300
3
2.5
2
£ 1.5
o
£ 1
0.5
0
-0.5
M
2 1
/ \ / V \ 1 \ T x
\ i u b lY \ V V.
100 200 L (mm)
300
Fig. 5. Distribution of flow pressure (a) and velocity (b) when opening the valve by: 1 — 1 mm, 2 14 mm.
2>
Fig. 6. Gas flow contours. _RUSSIAN JOURNAL OF NONLINEAR DYNAMICS, 2021, 17(1), 119 138
flow. With a disk lifting height of 1 mm, one can observe pressure and velocity jumps in the gap between the seat and the disk corresponding to the jet character of the gas flow with shock waves in this region (Fig. 6). The change in speed beyond the critical section is uneven in both cases.
To assess the influence of nonclassical flow regions on the device operation, we have performed calculations of the gas flow in a safety valve with the following coefficients of the Van der Waals equation of state: a = 0, b = 0, i.e., for ideal gas. We have calculated the ideality index Z, which is the ratio of the equations of state for ideal and real gases:
z = l
P
8 (p + 3p2)(l - p/3) p
Our comparison of the results showed that the presence of regions with negative derivative values affects the features of the gas flow, the main of which is the behavior of the local sound speed and the Mach number. Figure 7 shows the change in the ideality index Z and the Mach number along the axis of the cylindrical part and the disk surface generatrix. The valve opening is 14 mm. The greatest deviation from ideality is observed above the region of passing the sound speed.
N
4 t
3 7 y 1 ' PxC / W V *
I
2 " - ~ ~ l"----
100
200
L (mm)
Fig. 7. Distribution of the Mach number and the ideality index in the flow, 1, 2 gases, 3, 4 — Z for real and ideal gases, respectively.
300
M for real and ideal
A comparison of absolute pressures is shown in Fig. 8. The general flow pattern is preserved, but the difference in pressures on the disk surface leads to a difference in the total gas-dynamic forces. For an ideal equation of state, the force is 3,976 N, for a real one, 3,692 N. An 8% difference is significant when setting the valve for stable operation. The values of substance density in the flow differ more strongly (Fig. 9). The maximum difference is approximately 400 kg/m3 near the critical section. The most significant deviation is recorded in the calculated temperature values (Fig. 10). Along the entire length of the channel, the temperature difference between ideal and real media is more than 300 K.
Significant differences in the gas-dynamic characteristics of the flow for models of ideal and real gases confirm the need for making calculations using the real equation of state of the working substance. This allows us to correctly describe the dynamic processes in a safety valve at the stages of opening and closing, when an accurate determination of the gas-dynamic force of the disk is especially important.
20
15
e
10
1
2 \ \ _ /
V, r-------
100 200 L (mm)
Fig. 8. Pressure distribution, 1 — real gas, 2 — ideal gas.
300
1200
Ö u
P
900
600
300
2
1 ^ \ 1 V -^n i
— I"
V №—
100 200 L (mm)
Fig. 9. Density distribution, 1 — real gas, 2 — ideal gas.
300
800
~ 600 I
| 400
I
fS
M 200
1 ^______ -i^-V. 1-,---
2
100 200 L (mm)
300
Fig. 10. Temperature distribution, 1 — real gas, 2 — ideal gas.
The movement dynamics of the locking element (disk) is determined by the gas-dynamic force acting on the disk surface. The disk movement in the axial direction £, coinciding with the symmetry axis of the body's internal volume, is determined by the action of the gas force Ff, the spring elastic force Fs = Ks(£ + £0), the damping friction force Fd = Kdn, and the gravity msg. Here, Ks is the spring stiffness coefficient, £ is the disk travel (lifting height), £0 is the initial spring compression, Kd is the damping (friction) coefficient, and n is the disk travel speed. The disk travel equations are as follows:
dn
= Ff - Fs - Fd - msg,
d£
Tt=r]-
(5.3)
The problem of the joint solution to the gas-dynamic equations (5.1), (5.2) and the disk travel equations (5.3) is conjugate. The connection is carried out through the gas-dynamic force, determined through the integral of pressure and viscous friction stress tw on the disk surface:
Ff = J(p + tw )rdrdd.
Nonstationary gas-dynamic phenomena are observed when the disk moves during the opening and closing of its flow section. There are possible modes of monotonic and oscillatory movement of the disk [25]. With full opening of the flow area, there are also nonmonotonic flow modes determined by the device design and flow conditions. Figure 4 shows vortex flow areas. These vortices are unsteady; they are associated with the propagation of sound waves and affect the type of disk motion.
The time variation of the gas-dynamic force acting on the disk, and the pressure at local points of the inlet area are shown in Fig. 11 and Fig. 12.
Figure 11 shows the fluctuations in the gas-dynamic force with a fully open flow section. The frequency of these fluctuations is approximately 100 Hertz. Our studies have shown that the inlet pipe of a safety device is also a source of pressure fluctuations when the substance flows out of the controlled container. Figure 12 shows the change in pressure during one cycle at different distances from the inlet to the device pipe. In the region upstream of the minimum
2600
2590
u u
o №
2580
2570
Fig. 11. Fluctuations in gas dynamic force with full valve opening.
1.05
a
U 0.9
W3 W
u
Jh fc
0.75
20 25 30 35
t (ms)
Fig. 12. Pressure fluctuations at different distances from the inlet, 1 — 1 mm, 2 — 80 mm, 3 — 160 mm.
flow section, waves are formed with a frequency of about 50 Hertz. In the same area, one can see a transition from a nonclassical to a classical flow regime associated with moving away from the critical isotherm for the real state of substance.
6. Conclusion
In this study, we have built a mathematical simulation of a real inviscid gas flow using the Van der Waals equation of state. Based on the Godunov method, we have obtained an algorithm for the digital calculation of real gas flows with a local approximation of the Van der Waals equation of state by a two-term equation of state.
Our test calculations for solving Riemann's discontinuity disintegration problem have shown that the resulting algorithm is applicable both in the regions of classical gas behavior and in those with nonclassical behavior at negative fundamental derivative values. We have demonstrated that the presence of regions with nonclassical behavior of gas leads to blurring of the shock wave front.
Based on the developed Godunov scheme for solving gas-dynamic equations in a three-dimensional setting, we have carried out a digital simulation of the safety valve operation for systems filled with real gas at a disk lift height of 1 mm and 14 mm. We have determined the regions of classical and nonclassical behavior of gas, simulated the gas flow in a safety valve for ideal and real equations of state. We have shown that a correct description of the dynamic processes in the device at the stages of its opening and closing is possible only on the basis of a real equation of state.
The joint solution to the equations of real gas dynamics and the equations for the shut-off element travel showed that the disk travel during the working medium discharge generates unsteady gas-dynamic phenomena in the flow, whereas the device inlet pipe is a source of pressure fluctuations.
References
[1] Argrow, B.M., Computational Analysis of Dense Gas Shock Tube Flow, Shock Waves, 1996, vol.6, pp.241-248.
[2] Beune, A., Analysis of High-Pressure Safety Valves, PhD Thesis, Eindhoven: Technische Universiteit Eindhoven, 2009.
[3] Casari, N., Pinelli, M., Suman, A., and Pavanelli, G., Reducing Pressure Valve with Real Gases: An Integrated Approach for the Design, in 73rd Conf. of the Italian Thermal Machines Engineering Association, 2018, pp. 607-614.
[4] Colella, P. and Glaz, H. M., Efficient Solution Algorithms for the Riemann Problem for Real Gases, J. Comput. Phys, 1985, vol. 59, no. 2, pp. 264-289.
[5] Colonna, P. and Guardone, A., Molecular Interpretation of Nonclassical Gasdynamics of Dense Vapors under the vanderWaals Model, Phys. Fluids, 2006, vol. 18, no. 5, 056101, 14 pp.
[6] Dahmen, W., Muller, S., and Vob, A., Riemann Problem for the Euler Equation with Non-Convex Equation of State including Phase Transitions, in Analysis and Numerics for Conservation Laws, G. Warnecke (Ed.), Berlin: Springer, 2005, pp. 137-162.
[7] Davidson, R., Pre-Combustion Capture of CO2 in IGCC Plants, London: IEA Clean Coal Centre, 2011.
[8] Fossati, M. and Quartapelle, L., The Riemann Problem for Hyperbolic Equations under a Nonconvex Flux with Two Inflection Points, arXiv:1402.5906 (2014).
[9] Numerical Solution of Multidimensional Problems of Gas Dynamics, S.K.Godunov (Ed.), Moscow: Nauka, 1976 (Russian).
[10] Guardone, A. and Vigevano, L., Roe Linearization for the vanderWaals Gas, J. Comput. Phys., 2002, vol. 175, no. 1, pp. 50-78.
[11] Guardone, A. and Vimercati, D., Exact Solutions to Non-Classical Steady Nozzle Flows of Bethe-Zel'dovich-Thompson Fluids, J. Fluid Mech, 2016, vol.800, pp. 278-306.
[12] Henderson, L. F., General Laws for Propagation of Shock Waves through Matter, in Handbook of Shock Waves: Vol.1. Theoretical, Experimental, and Numerical Techniques, G.Ben-Dor, O.Igra, R. Elperin (Eds.), San Diego: Acad. Press, 2000, pp. 143-183.
[13] Kopyshev, V. P., Medvedev, A. B., and Khrustalev, V. V., Equation of State of Explosion Products on the Basis of a Modified vanderWaals Model, Combust. Explos. Shock Waves, 2006, vol. 42, no. 1, pp. 76-87.
[14] Kulikovskii, A. G., Pogorelov, N.V., and Semenov, A. Yu., Mathematical Problems of the Numerical Solution of Hyperbolic Systems of Equations, Chapman & Hall/CRC Monogr. Surv. Pure Appl. Math., vol. 118, New York: Chapman & Hall/CRC, 2001.
[15] Kuzenov, V.V., Ryzhkov, S.V., and Starostin, A.V., Development of a Mathematical Model and the Numerical Solution Method in a Combined Impact Scheme for MIF Target, Russian J. Nonlinear Dyn., 2020, vol. 16, no. 2, pp. 325-341.
[16] 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: Part1, Vestnik YuUrGU. Ser. Mat. Model. Progr., 2019, vol.12, no. 1, pp. 32-43 (Russian).
[17] Lopez-Echeverry, J., Reif-Acherman, S., and Araujo-Lopez, E., Peng-Robinson Equation of State: 40 Years through Cubics, Fluid Ph. Equilibria, 2017, vol.447, pp. 39-71.
[18] Moiseev, N.Ya. and Muhamadieva, T.A., Newton's Method As Applied to the Riemann Problem for Media with General Equations of State, Comput. Math. Math. Phys., 2008, vol.48, no. 6, pp. 1039-1047; see also: Zh. Vychisl. Mat. Mat. Fiz., 2008, vol.48, no. 6, pp. 1102-1110.
[19] Osher, S. and Solomon, F., Upwind Difference Schemes for Hyperbolic Systems of Conservation Laws, Math. Comput, 1982, vol. 38, no. 158, pp. 339-374.
[20] Papes, I., Abdelli, L., Degroote, J., and Vierendeels, J., 3D CFD Analysis of a Twin Screw Expander With Different Real Gas Models for R245fa, in ASME 2015 International Mechanical Engineering Congress and Exposition (Houston, Tex., Nov 13-19, 2015), Paper No: IMECE2015-53388, V07AT09A045, 8 pp.
[21] Polezhaev, V.I., Methods for Modeling Convective and Wave Processes and Heat Transfer in Real-Critical Media: An Overview, Fluid Dyn., 2011, vol.46, no. 1, pp. 1-15.
[22] Prokopov, G.P. and Severin, A.V., Rational Realization of Godunov's Method, Preprint No. 29, Moscow: KIAM, 2009 (Russian).
[23] Quartapelle, L., Castelletti, L., Guardone, A., and Quaranta, G., Solution of the Riemann Problem of Classical Gasdynamics, J. Comput. Phys, 2003, vol. 190, no. 1, pp. 118-140.
[24] Raeder, T., Tenenev, V. A., and Chernova, A. A., Determination of Flow Characteristics in Technological Processes with Controlled Pressure, Instruments and Methods of Measurement, 2020, vol. 11, no.3, pp. 204-211.
[25] Raeder, T., Tenenev, V. A., and Chernova, A. A., Numerical Simulation of Unstable Operating Modes of a Safety Valve, Vestn. Tomsk. Univ. Mat. Mekh, 2020, no. 68, pp. 141-157 (Russian).
[26] Raeder, T., Tenenev, V., Chernova, A., and Koroleva, M., Multilevel Simulation of Direct Operarted Safety Valve, in 2018 Ivannikov Ispras Open Conference (ISPRAS, Moscow, 2018), pp. 109-115.
[27] Raeder, T., Tenenev, V., and Koroleva, M., Numerical Simulation of the Working Process in a Safety Valve with Additional Gas-Dynamic Coupling, Intellekt. Sist. Proizv., 2020, vol. 18, no. 3, pp. 118-126 (Russian).
[28] Reid, R. C., Prausnitz, J.M., and Sherwood, Th. K., The Properties of Gases and Liquids, New York: McGraw-Hill, 1977.
[29] Rinaldi, E., Pecnik, R., and Colonna, P., Exact Jacobians for Implicit Navier - Stokes Simulations of Equilibrium Real Gas Flows, J. Comput. Phys, 2014, vol.270, pp.459-477.
[30] Sexton, A., Modeling Real Gases and Liquids Using a Modified vanderWaals Equation of State, Electronic Theses and Dissertations, Paper 1301, 2004.
[31] Toroa, E., Castrob, C., and Leec, B., A Novel Numerical Flux for the 3D Euler Equations with General Equation of State, J. Comput. Phys., 2015, vol. 303, pp. 80-94.
