INCORPORATION OF FLUID COMPRESSIBILITY INTO THE CALCULATION OF THE STATIONARY MODE OF OPERATION OF A HYDRAULIC DEVICE AT HIGH FLUID PRESSURES Текст научной статьи по специальности «Физика»

Russian Journal of Nonlinear Dynamics, 2021, vol. 17, no. 2, pp. 195-209. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd210205

NONLINEAR ENGINEERING AND ROBOTICS

MSC 2010: 35Q30, 76-10, 76N15

Incorporation of Fluid Compressibility into the Calculation of the Stationary Mode of Operation of a Hydraulic Device at High Fluid Pressures

This paper is concerned with assessing the correctness of applying various mathematical models for the calculation of the hydroshock phenomena in technical devices for modes close to critical parameters of the fluid. We study the applicability limits of the equation of state for an incompressible fluid (the assumption of constancy of the medium density) to the simulation of processes of the safety valve operation for high values of pressures in the valve. We present a scheme for adapting the numerical method of S. K. Godunov for calculation of flows of incompressible fluids. A generalization of the method for the Mie-Gruneisen equation of state is made using an algorithm of local approximation. A detailed validation and verification of the developed numerical method is provided, and relevant schemes and algorithms are given. Modeling of the hydroshock phenomenon under the valve actuation within the incompressible fluid model is carried out by the openFoam software. The comparison of the results for the weakly compressible

Received April 22, 2021 Accepted June 13, 2021

This work was performed as part of 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 Sciences" AAAA-A19-119092690104-4.

Thomas Raeder [email protected] Alena A. Chernova [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. Tatiany Baramzinoi 34, Izhevsk, 426067 Russia

T. Raeder, V. A. Tenenev, A. A. Chernova

and incompressible fluid models allows an estimation of the applicability ranges for the proposed schemes and algorithms. It is shown that the problem of the hydroshock phenomenon is correctly solved using the model of an incompressible fluid for the modes characterized by pressure ratios of no more than 1000 at the boundary of media discontinuity. For all pressure ratios exceeding 1000, it is necessary to apply the proposed weakly compressible fluid approach along with the Mie-Grüneisen equation of state.

Keywords: hydraulic device, mathematical model, numerical simulation, Godunov's method, Mie - Grüneisen equation of state, water, weakly compressible fluid approach, incompressible fluid

1. Introduction

Stop valves of various configurations which are, in fact, hydraulic devices of various complexity are significantly widespread in industry, housing and public utilities, and gas and oil industry. Working processes in such devices [1] are associated with the hydroshock phenomena and are characterized by high values of pressure. The most common working substance for this kind of devices is water with its compressibility for modes close to critical [3], causing incorrect results of applying the common incompressible fluid approach [17, 23-25].

Refs. [2, 20, 22] present numerical studies of nonstationary working processes in hydraulic devices with a compressible gas used as a working substance. The studies presented in these papers assume the application of the Euler and Navier-Stokes equations for a compressible medium using the Mendeleev-Clapeyron equation of state for a perfect gas within the method of S. K. Godunov and show good agreement of the calculated data with the experimental ones. Thus, in [2, 20] the original algorithm of the solution of nonstationary spatial equations of gas dynamics for a perfect gas is proposed. This approach requires substantial adaptation to be used to study the dynamics of a compressible fluid. For example, the transition from a perfect gas to a viscous fluid requires correcting the mathematical model, in particular, the transition from the Mendeleev - Clyperon equation of state to rheological equations describing the behavior of the medium (water) as the critical state is approached.

The standard model of an incompressible fluid defined by some value of the constant density of the substance does not take into account the compressible effects, which, at high operating pressures in hydraulic devices, may lead to a discrepancy between the calculated and experimental data. In addition, a significant increase in the medium pressure in the hydraulic device before its actuation leads in some cases [20] to the transition of parameters of the working media to the range of critical values. In this regard, when modeling such modes of hydraulic device operation, it is necessary to apply models which are correct for states of the substance close to critical

In [5-7], the problems of applying the Van der Waals equation of state of a real gas and its modifications (Redlich - Kwong and Peng-Robinson equations) are addressed to describe the medium behavior under conditions close to critical ones. The best known equation of state for a real gas is the Van der Waals equation, which takes into account intermolecular interactions and the real (nonzero) intrinsic volume of molecules. This model is able to reflect a real set of different gas/liquid states depending on the distance from the critical point. Thus, on the basis of the Van der Waals polytropic gas model, Ref. [21] studies the nonclassical gas-dynamic behavior of vapors near the liquid-vapor saturation curve. A modification of the Van der Waals equation is used to describe explosives in the region of pressures and temperatures typical for detonation processes [7]. In [4] it is shown that the Jones-Wilkins-Lee equation is used to

study explosive processes at high pressures, and the equation in the form of Mie- Grüneisen is applicable to analyze the fluid flow at high pressures. In accordance with [8], the Jones-Wilkins-Lee equation of state finds wide application in analyses of various explosive and, in particular, detonation processes.

Since the functioning of hydraulic devices, including direct operated safety valves, involves the implementation of large pressure drops in the systems under study, the Mie-Grüneisen equation of state is of particular interest for describing the rheology of a weakly compressible fluid for critical states.

In [9], a numerical algorithm is presented for the calculation of intense shock waves in multiphase condensed medium problems, based on the joint use of the conservative method of S.K.Godunov for high-order accuracy in the Eulerian form [10] and a thermodynamic model based on the Mie-Grüneisen equation of state along with the Hugonio equation. This thermo-dynamic model is applied to a wide range of nonporous condensed phases where the multiphase process is modeled by designing a single "effective phase" whose density, specific energy and elastic properties are determined by the average values of individual phases. To decompose the "effective phase" flows back into their respective individual components, a second-order accuracy algorithm for reconstructing the interphase boundary is used. A combination of a two-dimensional version of the splitting operator and an adaptive mesh refinement algorithm is used in [9] to solve the shock-wave geophysics problems.

Another approach to the application of the equations of state of real media in solving the problem of breakdown of an arbitrary discontinuity by the Newton method is described in [11]. The algebraic isentropic equation and some explicit solutions are written out in explicit form for media governed by the Mie-Grüneisen equation of state. An algorithm is proposed for the approximate calculation of Riemann integrals that has a higher accuracy and speed of convergence than the algorithm in the Simpson method. Possibilities of the approach are illustrated by examples of solving problems of an arbitrary discontinuity breakdown for media governed by the Mie-Grüneisen equation of state.

The paper [12] summarizes the current state of thermophysical properties of pseudo-pure and pure fluids. More accurate thermodynamic properties of fluids based on multiparameter explicit-type formulations taking into account the Helmholtz energy have been obtained. A larger set of methods has been used to estimate the transfer properties, including an extended method for corresponding states. This method can be used for fluids that do not have specific relationships for the transfer properties.

In [13], a modification of the flow splitting method proposed in [12] is considered. The original scheme was developed for one-dimensional Euler equations within an ideal gas model. The proposed modification includes three significant differences from the original scheme: 1) solving three-dimensional Euler equations on arbitrary meshes; 2) using the general equation of state; and 3) using the semidiscrete ADER (Arbitrary Accuracy Derivative Riemann) problem method on an arbitrary mesh to improve the order of accuracy for numerical schemes in time and space. The finite method is evaluated for accuracy, stability, and efficiency for a number of test problems. The accuracy is evaluated by investigating the convergence rate of schemes up to the 4th degree of accuracy in time and space on unstructured meshes.

Thus, this paper addresses both problems of extending the applicability of the method of S. K. Godunov for low compressible viscous fluids and assessment of the correctness of application for various mathematical models to the analysis of the hydroshock phenomena in technical devices for modes close to critical parameters of the fluid.

2. Mathematical Models

We consider the flow of water in a hydraulic device (safety valve) with the spatial discretization of internal channels shown in Fig. 1.

Fig. 1. Discretization over the space of the inner cavity of the hydraulic device

A physical concept of instantaneous change of working parameters of the fluid passing through a narrow gap region (area A, Fig. 1) is used. The solution of the problem of an arbitrary discontinuity breakdown is developed using two approaches: taking into account and without taking into account the fluid compressibility. Mathematical models, numerical schemes and algorithms used within the chosen approaches are presented in detail in Table 1.

3. Testing of calculation schemes and algorithms

To test the proposed numerical method for modeling of weakly compressible flows at high pressure drops, the problem of an arbitrary discontinuity breakdown [20] with an explicit solution is considered.

The Riemann problem involves solving the system of one-dimensional Euler gas-dynamic equations in the form

d G d H

— + — = 0, 3.1

dt dx

where t is time and x is a coordinate. Vectors G and H are defined as

G

P

pu e

H

pu

2

pu + p (e + p) u

where p is the medium density, u is the velocity, e = p(e + u2/2) is the total energy per unit volume of gas, and e is the internal energy per unit mass assigned by the equation of state.

Table 1. Characteristics of the models, schemes and algorithms used

Model of the weakly compressible fluid (water) Model of the incompressible fluid (water)

Mathematical model

t + Vpu = 0, p^ = pF- Vp + divP, p jjE = pFu + V(Pu) + Vq, Vu = 0, ^ + V • (uu) =g-Vp + V- (¡/Vu) , p = const,

(7-1 )P Pk - (du' | dxj I dxj dxj, I ' 7 = 1.638; p = 5; pk = 1180.5 kg/m3; cfc = 2090.1 m/s; Cv = 3279 Dj/kgK for p e [10; 1000] MPa; Ta G [293; 448] K; pe [1004; 1180] kg/m3

Initial conditions

up to the narrow gap: p = 0.1319, p = 1, u = v = 0; behind the gap: p = 0.1319-10"4, p = 0.9, u = v = 0. up to the narrow gap: p = 1 GPa, u = v = 0; behind the gap: p = 0.1 MPa, u = v = 0.

Boundary conditions

— nonflow condition on solid surfaces; — symmetry condition on the axis; — parameter extrapolation at the outlet boundary; — enthalpy of braking and constant entropy condition are assigned at the inlet boundary [1]. — nonflow condition on solid surfaces; — symmetry condition on the axis; — soft boundary conditions at the outlet boundary; — total pressure is assigned at the inlet boundary.

Numerical methods, schemes and algorithms

The solution of the resulting system of equations is framed by S.Iv. Godunov's method [10, 14]. In this case the local approximation of the equation of state allows one to find the initial approximation of variables for solving the resulting system of equations. The solution of the Navier - Stokes equations is framed by the method of control volumes using the combined algorithm Pimple in the open package OpenFOAM, solver pimpleFoam.

The Riemann problem involves solving equations 3.1 under the initial conditions

^. . (G\,x < 0, .

G(0'x) = {G2,x< 0, (3-2)

In this case, we speak of the discontinuity breakdown problem [14]. The structure of the flow is determined by a combination of rarefaction waves and shock waves with a contact discontinuity or vacuum area. The characteristics of the flow are fully described by the known relations for the shock wave and rarefaction wave. For the shock wave, the jump of parameters p, p, u, e is

determined by the Rankin -Hugonio conditions [6]:

[p]

M =

[p]

1

p

[e] = -0.5(P + p)

(3.3)

(3.4)

(3.5)

where

= U - u, [p] = P - p, [e] = E - e,

= — i, u, p, e, p are the velocity pressure,

internal energy, and density, respectively, before the shock wave front, and U, P, E, R are the velocity, pressure, internal energy, and density behind the shock wave front; a is the mass velocity.

For a rarefaction wave, there is an isentropic equation in the form

dp 2 dp

and an equation corresponding to the constancy of the Riemann invariants

du 1 dp pc

(3.6)

(3.7)

where c is the speed of sound. For the rarefaction wave, the Poisson adiabatic equation is also written [5]

^ = (3-8)

dp p2

In addition, the continuity condition at the contact discontinuity of velocity and pressure is used.

Taking into account the conditions considered in [3, 4, 20], we compose a system of four equations to determine the values of P, U, R\, R2 (where R\ and R2 are the densities of the medium before and after the discontinuity boundary, respectively). The transformation of relations (3.3) followed by the exclusion of U from them leads to the form

P2

pi

Wi W2

- ui + u2 = 0,

(3.9)

where the expressions for the mass velocity derived from equations (3.4) and (3.7) are determined by the expression

' D \/P ~ Pi D ^ RiPi~5-—> P>Pi,

Ri - pi Wi ={ P - pi

Pi

f

J pc

P < pi

(3.10)

and correspond to (3.3) for the shock wave P > pi or the Poisson adiabat (3.8) for rarefaction waves.

For shock waves, the energy equation is

E(P,Ri) + 0.5(P+Pi) ( -I-i ) = 0,

Ri pi,

(3.11)

a

where the function E(P, R) corresponds to the equation of state of the matter considered. For rarefaction waves, it is necessary to numerically integrate the differential equation (3.8) if the rheology of the studied medium is described using the complex equation of state. Then, using the Mie-Gruneisen equation [11], the function e(p,p) in (3.8) can be defined as

p - k4

p — i

e =

(7-!)(*-!) + V^1-«")

(Y — 1)p

(3.12)

where 7, T, d, A, B, a, ft, pk, Ck are some constants defined for these types of equations, 5 = -j^,

and equation (3.6) is integrated in quadratures.

Then the expressions for the pressure and sound velocity changes in the rarefaction wave for the Mie-Gruneisen equation [13] are

PslPi, Pi, p) = - pk( 1 + f - P), Pk = l''/i

Y

a(S ) = Y

Pi + Pk (1 + — 6»)

Pi

c = Y

P + Pk [1 + (P — Y )6»

P

(3.13)

(3.14)

(3.15)

Let us write the second and third equations for the unknowns G as

/2(G) = 0 =

/3(G) = 0 =

E(P,Ri ) + 0.5(P + pi)

Ri

P — ps(pi,pi,Ri), P > Pi.

E(P,R2 ) + 0.5(P + P2)

R2

P>Pi,

P>P2.

P - Ps(P2 ,P2,R2), P > P2-The value of the velocity is determined by the condition

P - Pi

(3.16)

(3.17)

U = Ui±

Wi

The system of nonlinear equations

fj (G) = 0, j = 1,2,3

(3.18)

(3.19)

is solved at some initial approximation Go by the Newton method. When calculating the mass

p

velocity in the case of a rarefaction wave, it is necessary to calculate the integral J In quadra-

p- p

tures, this integral can be written only for some particular cases [11].

The explicit solution of the Riemann problem in the case of equations of state for a real gas is associated with an iterative solution of the system of equations (3.19), which is an expensive, in terms of computational costs, operation within the numerical implementation difference schemes of Godunov type. The approach of local approximation of the equation of state by a binomial

P

1

1

equation [11, 15, 16] is known. If the equation of state is assigned as 0(e, p) = 0, it can be repre-

sented in the form of a binomial equation 0(e, p) = e — determined locally based on conditions [16]:

P + TPo- Pc% Pi 7-1)

with parameters Y, po, ck,

dp

, , 9</> P + 7Po-PTp

(£'P)=0; (7 — l)p2 =0;

de

= 1 -

dp

lie

(y -1)

0.

(3.20)

The local approximation of parameters 7, p0, ck for the Mie - Gruneisen equation can be represented as

Po = Pk\ 1 +

At - 7 ( P_

(3.21)

Y \Pk

For numerical implementation within the author's code it is convenient to represent the equations of motion in dimensionless form by means of the ratios: flow velocity to the critical

speed of sound ü = density to the critical density 0 = pressure in the form p = c8

p .

9 >

p*c4

internal energy in the form e =

For the Mie-Grüneisen equation the density is related to the value obtained from the solution of the equation

poo

e(Too) =

2

PkC% p — 1

(7-D(f-l) +

Y — P

i-m1

(Y — 1)Poo

(3.22)

where p00 and Too are the known braking parameters for the flow.

P00 + Pk

p00

The pressure is related to the value of PooCqo> where coo = \ 7-

the internal energy to c^o, and the velocity to coo.

Let us consider the possibility of using the Mie-Grüneisen equation using the water flow as an example. In [17], equations describing the dependences of the volume on pressure and temperature are given as a polynomial and an exponent for pressures up to 1000 bar. Tabulated dependences also up to 1000 bar are contained in the reference book [18]. A more complex wide-range equation of state for water and vapor in analytical form is developed in [19]. Based on these data, parameters of the Mie-Grüneisen equation were chosen by the least-squares method Y = 1.638; p = 5; pk = 1180.5 kg/m3; ck = 2090.1 m/s; Cv = 3279 dj/kgK, valid for the ranges of p e [10; 1000] MPa; T e [293; 448] K; p e [1004; 1180] kg/m3.

The approximation diagram of the dependence p(p, T) is shown in Figure 2.

Testing of the proposed computational algorithm involves comparing the numerical solution of the Riemann problem in a shock tube with the explicit solution. Three design cases of the problem of an arbitrary discontinuity breakdown are considered: propagation of two rarefaction waves, propagation of a rarefaction wave to the left and a shock wave to the right simultaneously, and propagation of two shock waves. The initial parameters in the shock tube are given in Table 2.

A comparison of the explicit solution with the numerical one is shown in Figure 3.

Figure 3a shows that the results of modeling the breakdown of an arbitrary discontinuity within the framework of a real weakly compressible fluid model (the Mie-Grüneisen equation

ß

P

Fig. 2. Pressure vs. density and temperature determined by the Mie-Gruneisen equation

Table 2. Initial parameters in the shock tube

Pi U1 Pi P2 Uo Pi

two rarefaction waves

0.13 -0.1 1 0.13 0.1 0.9

rarefaction wave on the left, shock wave on the right

0.13 0 1 0.0013 0 0.9

two shock waves

0.13 0.2 1 0.0013 -0.2 0.9

of state) are in good qualitative agreement with the results of the explicit solution in the areas before and after the discontinuity. For all curves, two rarefaction waves are registered, however, in the area of the flow discontinuity, blurring of the rarefaction wave front is observed. In this case, the analysis of test calculations for flow configurations with a shock wave and a rarefaction wave (Fig. 3b) and two shock waves (Fig. 3c) shows that only the rarefaction wave front is subject to blurring, therefore, the registered differences between the numerical and explicit solutions are caused not by the influence of the scheme viscosity, but by the transition of fluid parameters from the critical area to the classical one. A quantitative comparison of the calculated parameters P, U, Ri, R2 is presented in Table 3.

It follows from the presented test calculations that the local approximation principle is applicable for calculation of weakly compressible fluids by S. K. Godunov's method along with the Mie-Gruneisen equation of state. The test calculations confirm the correctness and adequacy of the proposed algorithm for calculating the flows of weakly compressible fluids by S. K. Godunov's method.

1.2 1 0.8 0.6 0.4 0.2 0

— U-exac — ro-exac p-exac » U-calc ° ro-calc ° p-calc

-0.2 -0.6

-0.4 -0.2

0

x

(a)

1.2

1

s

* s 0.8

* Ci 0.6

va.

(M 0.4

«

Q. a, 0.2

H

0

-0.2

— U-exac — ro-exac p-exac • U-calc » ro-calc ■ p-calc

0.2 0.4 0.6

-0.6 -0.4 -0.2

0

x

(b)

0.2 0.4 0.6

1.2 1 0.8 0.6 0.4 0.2 0

-0.2

U-exac — ro-exac p-exac U-calc ° ro-calc ° p-calc

-0.6 -0.4 -0.2

0

x

(C)

0.2 0.4 0.6

Fig. 3. Comparison of the explicit and numerical solutions for: a) two rarefaction waves; b) rarefaction wave and shock wave; c) two shock waves

Table 3. Parameters P, U, Ri, R2

P U Ri i?2

two rarefaction waves

Explicit 0.05241 -0.01167 0.9082 0.7855

Godunov 0.05251 -0.0101 0.9069 0.7905

rarefaction wave on the left, shock wave on the right

Explicit 0.06229 0.0757 0.9218 0.9831

Godunov 0.06243 0.0774 0.9204 0.9849

two shock waves

Explicit 0.2876 0.0678 1.1248 1.162

Godunov 0.2874 0.0687 1.1246 1.164

4. Modeling of working processes in the hydraulic device

We consider the possibility of using models of weakly compressible and incompressible fluids for the numerical study of the working processes occurring in the safety valve (Fig. 1) when it operates in various modes including critical ones.

Three design modes of the device operation are considered: at the pressure difference in the valve (Pin/pout) equal to 100 (braking pressure at the inlet to the design area — 10 MPa), 1000 (braking pressure at the inlet to the design area — 100 MPa) and 10000 (braking pressure at the inlet to the design area — 1 GPa). A comparison of changes in the estimates of integral characteristics of the valve obtained within the weakly compressible and incompressible fluid approaches is presented in Fig. 4.

P* p*

(a) (b)

Fig. 4. Integral characteristics (hydraulic force (a) and fluid flow (b)) of the safety valve, obtained within the weakly compressible (curves 2) and incompressible (curves 1) fluid models

Divergences of the estimates of the relative hydrodynamic force (F/p) fixed for the pressure drops (Pin/pout) of 100 and 1000 are connected with an error of the Mie - Gruneisen equation of state when using it to describe the behavior of weakly compressible fluids under standard (noncritical) conditions, and deviations in values of flow characteristics for the pressure drop (Pin/pout) of 10000 are due to the low applicability of the incompressible fluid model for the description of compressible flows in the critical area.

The pressure fields obtained within weakly compressible and incompressible viscous fluid (water) models for pressure drops (Pin/pout) of 10000 are shown in Figure 5, and the spatial flow tubes are shown in Figure 6. In order to ensure the comparability of the results when using different equations of state of the fluid, a further comparison of approaches to the numerical modeling of working processes in the safety valve is given strictly for critical modes, namely, for the pressure drop (Pin/pout) equal to 10000.

Figure 5 shows that the pressure fields obtained within the weakly compressible and incompressible fluid approaches are qualitatively and quantitatively similar, but a number of fundamental distinctive features should be noted. In particular, the application of the incompressible fluid approach for modeling the fluid flow at high pressures leads to an incorrect estimation of the hydrodynamic characteristics in the area below the narrow gap; it causes the physically unreasonable jump-like type of the pressure field under the locking disk (Fig. 5b) and local dis-

0.0e+00

(a)

(b)

Fig. 5. Pressure fields in the hydraulic device obtained within the weakly compressible (a) and incompressible (b) fluid models

(b)

Fig. 6. Flow tubes in the safety valve, obtained within the weakly compressible (a) and incompressible (b) fluid models

turbances near the rod. However, configurations of the pressure fields in the compression zones (up to the narrow gap) and in the critical section for both approaches coincide.

As Fig. 6 shows, the structure of fluid flow in the valve in the critical mode (pressure drop (Pin/pout) of 10000, water pressure up to the narrow gap 1 GPa) calculated within the incompressible fluid model and that calculated within the weakly compressible fluid model are in good agreement with each other in the axisymmetric part of the valve. However, the fluid flows realizable in the outlet expanding channel of the valve are significantly different. The flow structure obtained in the incompressible fluid model is characterized by the presence of vortices of different configurations near the impermeable surfaces of the socket, and the realizable flow loses its

0.2-1

— wisCompr u

0.4 n

— wisCompr

— nonCompr

0.1

* Q.

0.05

-0.05

0.15

0

0

50 100 150 200 250 0

50 100 150 200 250

(a)

(b)

Fig. 7. Distribution of relative values of pressures (a) and velocities (b) along the rod surface for the angle 0°

stability after 1 ms of calculation. Also, it should be noted that there is some curvature of the flow tubes for the incompressible fluid model near the outlet line.

In order to clarify the limits of applicability of the mathematical models, an additional quantitative comparison of the calculated data has been performed. So, Fig. 7 shows the changes along the flow (at the inner boundary of the area in the plane 0°) of relative pressures and velocities obtained within the weakly compressible and incompressible fluid models.

It should be noted that the distribution of pressures obtained within the incompressible fluid approach (Fig. 7, red curve) is characterized by the presence of negative values after the flow passes a narrow gap. That is, the application of the incompressible fluid approach for modeling of flows with essential pressure drops does not allow calculation of correct parameters of the fluid behind the wave front. Fluctuations of pressure curves in the area behind the wave front are connected with topological features of the flow structure in this area. It should be noted that the pressure distributions obtained within the weakly compressible fluid approach using the Mie - Gruneisen equation of state (Fig. 7a, blue curve) contain no areas of negative pressures. The resulting velocity distributions (Fig. 7b) are qualitatively and quantitatively similar, but the application of the incompressible fluid approach for critical flow conditions causes underestimation of the flow velocities.

Figure 8 shows relative velocity profiles in a narrow gap for a realizable pressure drop (Pin/pout) of 10000 (the braking pressure at the entrance to the calculation region is 1 GPa).

It should be noted that the relative velocity profiles in the narrow gap for the pressure drop (Pin/pout) of 10000 are qualitatively and quantitatively similar. The asymmetry of the relative velocity profile can be noted within both approaches.

The combination of all above-mentioned factors indicates that the incompressible fluid model considered in this paper is not applicable to the calculation of hydrodynamics in modes close to critical ones, which are characterized by compressibility effects.

5. Conclusion

As a result of the work, an algorithm of adapting the Godunov method for calculation of weakly compressible fluids taking into account their real rheological properties within the Mie - Gruneisen equation of state has been proposed, tested and verified.

ln

0.9-

0.8-

0.70.6-

s

0.5-

0.40.30.20.10-

147

148

149

150

151

152

Fig. 8. Relative velocity profile in the narrow gap

Based on test calculations, it has been concluded that it is possible to apply the Godunov method with the local approximation of the equation of state for calculation of the fluid flow in safety valves under high pressures taking compressibility into account.

It is shown that it is expedient to perform modeling of working processes in hydraulic devices of the type considered within the incompressible fluid approach for calculations in standard modes (noncritical, with a pressure drop (Pin/pcut) up to 1000). To simulate the operation of the device in critical and transient modes (at pressure drops (Pin/pcut) of more than 1000 or in ultrahigh pressure modes), it is necessary to apply the proposed weakly compressible fluid approach involving the Mie-Gruneisen equation of state.

References

[1] Song, X., Cui, L., Cao, M., Cao, W., Park, Y., and Dempster, W. M., A CFD Analysis of the Dynamics of a Direct-Operated Safety Relief Valve Mounted on a Pressure Vessel, Energy Convers. Manag., 2014, vol.81, pp. 407-419.

[2] 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).

[3] Ismagilova, D. F., Ismagilova, R. F., and Tselischev, V. A., Mathematical Modeling of Water Hammer Protection System, Vestn. UGATU, 2014, vol.18, no.4(65), pp. 72-78 (Russian).

[4] Jin, Zh., Wei, L., Chen, L.L., Qian, J.Y., and Zhang, M., Numerical Simulation and Structure Improvement of Double Throttling in a High Parameter Pressure Reducing Valve, J. Zhejiang Univ. Sci. A, 2013, vol. 14, no. 2, pp. 137-146.

[5] 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.

[6] 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.

[7] Kopyshev, V. P., Medvedev, A. B., and Khrustalev, V. V., Equation of State of Explosion Products on the Basis of a Modified vander Waals Model, Combust. Explos. Shock Waves, 2006, vol. 42, no. 1, pp. 76-87.

[8] Trzcinski, W. A., Szymanczyk, L., and Kramarczyk, B., Determination of the Equation of State for the Detonation Products of Emulsion Explosives, Cent. Eur. J. Energ. Mater., 2019, vol. 16(1), pp. 49-64.

[9] Miller, G.H. and Puckett, E.G., A High-Order Godunov Method for Multiple Condensed Phases, J. Comput. Phys, 1996, vol. 128, no. 1, pp. 134-164.

[10] Godunov, S.K., A Difference Method for Numerical Calculation of Discontinuous Solutions of the Equations of Hydrodynamics, Mat. Sb. (N.S.), 1959, vol. 47(89), no. 3, pp. 271-306 (Russian).

[11] 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.

[12] Bell, I.H., Wronski, J., Quoilin, S., and Lemort, V., Pure and Pseudo-Pure Fluid Thermophysical Property Evaluation and the Open-Source Thermophysical Property Library CoolProp, Ind. Eng. Chem. Res., 2014, vol. 53, no. 6, pp. 2498-2508.

[13] Toro, E.F., Castro, C.E., and Lee, B.J., A Novel Numerical Flux for the 3D Euler Equations with General Equation of State, J. Comput. Phys., 2015, vol. 303, pp. 80-94.

[14] Numerical Solution of Multidimensional Problems of Gas Dynamics, S.K. Godunov (Ed.), Moscow: Nauka, 1976 (Russian).

[15] Prokopov, G.P. and Severin, A.V., Rational Realization of Godunov's Method, Preprint No. 29, Moscow: KIAM, 2009 (Russian).

[16] Kulikovskii, A. G., Pogorelov, N. V., and Semenov, A.Yu., Mathematical Aspects of Numerical Solution of Hyperbolic Systems, Monogr. Surv. Pure Appl. Math., vol. 118, Boca Raton, Fla.: Chapman & Hall/CRC, 2001.

[17] Likhachev, E.R., Equation of State of Liquid Mercury, Vestn. VGU. Ser. Fiz. Matem., 2014, no. 3, pp. 41-48 (Russian).

[18] Rivkin, S.L. and Alexandrov, A. A., The Thermophysical Properties of Water and Steam, Moscow: Energoatomizdat, 1984 (Russian).

[19] Nigmatulin, R. I. and Bolotnova, R. H., Wide-Range Equation of State for Water and Steam: Calculation Results, High Temp., 2008, vol. 46, no. 3, pp. 325-336; see also: Teplofiz. Vys. Temp., 2008, vol.46, no. 3, pp. 362-373.

[20] 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.

[21] 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.

[22] 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. 118126 (Russian).

[23] Borisov, A. V., Kilin, A. A., and Mamaev, I. S., The Dynamics of Vortex Rings: Leapfrogging, Choreographies and the Stability Problem, Regul. Chaotic Dyn., 2013, vol. 18, nos. 1-2, pp. 33-62.

[24] Kuzenov, V.V. and Ryzhkov, S.V., Mathematical Modeling of Plasma Dynamics for Processes in Capillary Discharges, Russian J. Nonlinear Dyn., 2019, vol. 15, no. 4, pp. 543-550.

[25] Mamaev, I. S., Tenenev, V. A., and Vetchanin, E. V., Dynamics of a Body with a Sharp Edge in a Viscous Fluid, Rus. J. Nonlin. Dyn., 2018, vol. 14, no. 4, pp. 473-494.