Simulation of gas-dynamic, pressure surges and adiabatic compression phenomena in geometrically complex respirator oxygen valves

Abstract A three-dimensional remeshed smoothed particle hydrodynamics method for the simulation of isotropic turbulence is used, the method is coupled with Brinkman penalisation technique for flow simulation inside the complex valve geometry. Simulations are carried out for three different pressure reserve quantities, to replicate the opening of the valve, two time based pressure inlet boundary condition functions were simulated along with an impulsively started scenario. A geometrical sensitivity analysis is provided, where the simulation is performed on a modified valve design which exhibits a damping effect on the gas dynamics and flow characteristics, which has a favourable effect on the valve functionality and safety. It is found that the capacity of the pressure reserve has a considerable effect on the simulated flow fields (velocity, temperature), as the temperature could rise 6.0 × the reference temperature, and up to 2.7 × the reference velocity. The numerical results are compared with a previous study carried out Rotarex S.A., demonstrating that the remeshed particle-mesh method coupled with Brinkman penalisation provides a good quality simulation and the results are in agreement with the reference solution.


Introduction
Industrial valves are widely used in a range of industrial applications like ultra high purity, medical, automotive, refrigerators, oxygen pipe lines, etc. Hence, the manufacturing companies require proper physical predictive models and smart efficient numerical methodologies to resolve complex flow problems including pressure surge, material erosion, adiabatic processes,.
Medical oxygen regulators/valves are required in respiratory therapy to improve the device's safety as the pressure regulator/valve reduces cylinder pressure and calibrates the gas inflow. In medical applications manufacturing, the main priorities are precision and safety. The precision of the regulator improves the gas delivery safety through better outlet pressure and stability of the inflow. In this paper we are presenting a numerical simulation of valve-opening gas dynamics in a medical oxygen valve. The study is conducted in collaboration with Rotarex S.A. which is a major medical gas valves and regulators manufacturer.
Due to cost cutting and global competition the market demands medical Oxygen systems with higher pressurised gas reserve, which are attached to the Oxygen valve, knowing that the current pressure in the cylinder is required to go as high as 300 bar and even above in the future. The difference between the pressure in the cylinder and the pressure in the valve can cause pressure surge and adiabatic compression phenomena as a result of the sudden change in gas velocity. Depending on the reserve pressure in the valve system, the materials which the valve is made, the occurrence of pressure shock wave and adiabatic compression ignition energy in the oxygen valve system can be triggered.
Gas-dynamic pressure surges and adiabatic compression are challenging issues when designing a butterfly valve. Depending on the field of application, the pressure surges can lead to considerable malfunctions, and in very critical cases to a total failure of the respective component as they are caused by an abrupt change in flow speed (as a result of the pressure difference) and can reach pressure valve's higher than the usual pressure conditions in normal circumstances [1].
Vu and Wang [2] presented a three dimensional model of Oxygen flow in a valve, and indicated vortex formation near the opening of the valve. Francis and Betts [3] modeled incompressible flow structure in a pressure relief valve, and later they measured the pressure distribution of a commercial valve disc [4]. Ahuja [5] performed a series of multiphase flow computational simulation of control valves with sub model for grid adaptation. Moncalvo et al. [6] discussed the effect of the degree of discretisation on the turbulent flow characteristics in two safety valves using the commercial software Ansys CFX-flo. Beune et al. [7] modeled a safety valve using Ansys CFX and simulated the compressible flow dynamics at pressure up to 3600 bar. Schmidt et al. [8] presented a detailed experimental and theoretical study on high pressure valves, including an optimised new design.
Dempster et al. [9] presented a two equation turbulence model for the prediction of gas flow rate in a refrigerator industrial valve. Jung et al. [10] showed that surge modeling is important to safeguard against breaches in the distribution system integrity.
We present a gas dynamic respirator oxygen valve simulation using an integrated hybrid remeshed smoothed particle hydrodynamics method. The computational approach introduced by Obeidat [11][12][13] is extended to the simulation of gas-dynamic pressure surges and adiabatic compression phenomena in complex geometries of a gas valve. To resolve the produced high Reynolds number turbulent structures, the particle-based numerical approach is adapted to large eddy simulation for turbulence modeling. To highlight the effect of the reserved pressure (pressure in the container attached to the valve) we present a numerical analysis of three different reserved pressure, at 300, 150 and 75 bar, also we show the effect of the speed of opening the valve by simulating two scenarios of a time-based valve opening along with the impulsively started one.

Pressure surge effect
A pressure surge is a transient sudden rise or fall of pressure in a fluid-conveying structure (valve, pipe,etc.). Opening the valve rapidly causes sudden change in gas velocity due to the pressure difference between the inside of the valve and the pressure reserve. The following related phenomena arise as a result of the pressure surge: • A sudden increase in the pressure gradient and pressure discontinuities can occur, causing pressure oscillations, • A wave moves through the fluid medium at greater velocity than the local speed of sound, a so-called shock wave, • A rapid surge of the velocity and density causes turbulent flow characteristics. • A rapid increase in temperature of the fluid, so called adiabatic or diabetic compression.
If the pressure surge is in excess of the rated capacity of a valve, it can cause ruptures in the structure causing the valve to be damaged. Additionally, the rapid increase of the temperature can cause the gas to ignite (depending on the nature of the gas). The turbulent flow causes the gas state to undergo irregular fluctuation, which results in cavitation in the valve body and in the case of any contamination (e.g. residuae alloy particles) in the valve, a spark can form, causing the gas to ignite.
Accurate simulations of flow, in which shock waves and the turbulent flow are presented and interact dynamically, are challenging because of the contradictory requirements of numerical methods used to simulate turbulence, which must minimise any numerical dissipation that would otherwise overwhelm the small scales, and shock-capturing schemes, which introduce numerical dissipation to stabilise the solution.
Shock waves are localised phenomena and in order to represent them accurately on a computational grid, most numerical schemes depend on numerical dissipation, which results in spreading the shock over a few grid points. The drawback of using shock-capturing schemes in smooth turbulent flow regions is, that the numerical dissipation overwhelms the physical dissipation, which is what turbulence modelling seek to avoid.

High pressure oxygen systems
Oxygen is one of the most common oxidising gas, oxygen is nonflammable, but it supports combustion, It reacts vigorously with combustible materials, and enhances a fire or explosion which will generate a large amount of energy in a short time. With high concentration of oxygen (100% in our study) both the risk of ignition and rate of combustion increases, at high pressures the ignition temperature decreases and combustion rate increases. Furthermore as we present in Section 7.1 increasing the pressure in the case of adiabatic compression at the inlet has a dramatic effect on the predicted temperature inside the valve. The increase of the temperature makes the risk of ignition more likely, as it decreases the amount of energy to be added to start a fire. With sufficient pressure and ignition energy, nearly all substances can be made to burn in pure oxygen.
The most critical ignition sources in a high pressure valve are: 1. Shock wave generation due to the pressure surge, as mentioned in Section 2. 2. Particle impact due to contamination. 3. Contamination with grease or oil. 4. Electric arc.
In the current study source 1, and 2 are more likely to occur. The valve body is often made of brass-materials (copper alloy) as they are cost efficient, ability to resist corrosion and naturally a soft metal, however the ignition temperature of brass at high pressure is low ≈ 1073 K and even lower for the polymeric materials used for sealing (e. g. Neoprene, Nylon or Polyethylene), refer to Figs. 1 and 2.
It is crucial to understand the risk of ignition within the valve by understanding gas dynamics, pressure surge and adiabatic compression and how they effect both the rising temperature of oxygen within the valve and the time the temperature takes to relax under the ignition threshold.

Governing equations
The compressible flow is governed by the balance equations with Brinkman penalisation as presented in [11,12]. The resulting conservation of mass is then, and the conservation of momentum, and finally the thermal energy equation, where D⋄ DT = ∂⋄ ∂t +(u⋅▽)(⋄ ) denotes the material derivative, u i is the velocity, p is the pressure, ρ is the density, τ ij is the shear stress, μ is the dynamic viscosity, δ ij is the Kronecker delta, τ sgs is the sub-grid stress tensor, which is zero in case of direct numerical simulation, k is the thermal conductivity, c v the specific heat at constant volume and the ( − ) denote to a filtered quantity using a discrete filter as in Eq. (13).
The shear stress τ ij is defined as, The sub-grid stress tensor τ sgs is defined as: SGS stresses interact between the resolved (grid) scales and the unresolved (subgrid) scales. In the present work we model the turbulent sub-grid stresses using the standard Smagorinsky model [14], which is an eddy viscosity model based on the Boussinesq approximation and for comprisable flow is defined as: which relates τ sgs to the eddy viscosity ν sgs and the resolved-scale strain rate tensor S ij . The eddy viscosity is defined as, where C s is a non dimensional constant for which values ranging from  0.1 to 1.0 and have been suggested in literature [15], Δ is the model The resolved-scale strain rate tensor S ij as, For the penalised term [12] u (oq)i is the solid velocity, ϕ is the porosity, η = αϕ is the normalised viscous permeability. Note that 0 < ϕ≪1, and 0 < η≪1. χ is the solid mask, and defined as, To close the system (Eqs. (1) and (2)), the following equation of state is used, with, where γ is the ratio of the specific heat capacities, and R is the gas constant.
In the presented work the flow is uniquely described by the Reynolds number Re = Uρ 0 L/μ, and the Mach number Ma = U/c. L is the characteristic length, ρ 0 is the reference density, U is the reference velocity, and c is the speed of sound.
The non-dimensional variables are obtained from the physical variables as, where the superscript (*) and the subscript (0) indicate the nondimensional and the reference quantities. The quantities L 0 , ρ 0 , U 0 , T 0 denote the characteristic length, density, velocity, temperature. The discrete filter applied to the quantities is chosen as a Gaussian filter of fourth order, where ∊ is the ratio of the mesh size to the cut-off length scales of the targeted filter and f i is the filtered quantity, for more details please refer to [16,13]

The hybrid remeshed smoothed particle method (hrSPH)
The presented algorithm is based on the previously developed hybrid remeshed smoothed particle hydrodynamics method (hrSPH) by Obeidat and Bordas [11]. For the simulation of flow around/inside complex geometry, Obeidat et al. [12] incorporate a Brinkman penalisations technique with the hrSPH. The computational domain is implicitly penalised by a mask function marking the regions where the solid geometry is located. A penalty term is added to all balance equations (implicit penalisation), however in this study we are explicitly penalising the flow by a split-algorithm technique.
The hrSPH method is divided into three parts: 1. Computing the rate of change (a) Particle-mesh interpolation of the temperature, mass and impulse of the particle.
where P(m) denotes the set of particles that contribute to mesh node x m , h the particle spacing, W the high order kernel, u p the three velocity components u, v, w, m p the mass of the particle, T p is the temperature of the particle, x m the position of mesh node m, and x p the position of particle p.
where u s is the velocity of the solid wall of the geometry, for no slip boundary condition u s = 0, for moving wall boundary condition u s can take any desired value. ΔT ( 2. Updating the particles This part takes place on the set of particles, where the interpolated rate of change in velocity, mass and temperature are used to update the velocity, mass, temperature and position of the particles, here showing the temperature update and particle advection.
3. Remeshing the particles: In case of distortion and particle clustering due to the turbulent nature of the flow, interpolate the strengths of the particles to the mesh using the M ′ 4 interpolation function, generate a new set of the particles, and then interpolate the strengths back to the new set of particles.
The M ′ 4 introduced by [17], interpolates the strength of the particles to the mesh so that the strengths are redistributed onto the surrounding mesh nodes. The redistribution of particle quantities is achieved using the 3rd order M ′ 4 kernel which in one dimension is expressed as: where |x| is the distance of the particle to the mesh, h is the spacing of the new uniform set of particles.

Numerical setup
The computational domain of the studied simplified oxygen valve is represented in Fig. 3  The simulation is performed with 100 particle in each coordinate, 27 × 2.6 GHz Dell C6420, Intel Xeon Gold 6132 are used to run the simulation, the simulation is completed after 3 h and the data is collected then.

Boundary conditions and initial conditions
At the inlet we specify a pressure boundary condition, and a bulk flow rate at Mach number Ma = 1 when the flow is fully developed, on the solid walls a no-slip boundary condition is prescribed. Initially the outlet pressure and initial pressure is set to 1 bar. The simulated gas is Oxygen with constant viscosity, and the initial temperature of the valve body is set to 290 K.

Results
We study the effect of the developing pressure state in the pressurised cylinder reserve on the fluid's dynamic and thermodynamic characteristics.
We are particularly interested in investigating the effect of increasing the pressure reserve on the pressure surge, and adiabatic compression development by simulating the three dimensional flow for three different pressure inlet conditions at the inlet and three different timebased valve opening functions.
We note here that the focus of the study is to investigate how the fluid is flowing from the inlet to the outlet along with the effect of the pressure surge on the fluid characteristics (velocity, temperature, etc.), and we are not interested in simulating a fully developed flow.
The time histories of mean velocity, temperature and Mach number are obtained from four points (locations) as illustrated in Fig. 3(a) and compared with a reference solution conducted by Rotarex S.A. using OpenFoam-sonicFoam solver and k − ω − SST turbulence model. In Section 7.1 we are presenting three cases with different pressure inlet boundary conditions at 75 bar, 150, and 300 bar. In this case the valve is open and the flow is impulsively started. In Section 7.2 we study the effect of opening the valve with respect to a time based opening functions, where the pressure at the inlet increases with time, mimicking a real valve opening. We characterise the flow with the following where ν is the kinematic viscosity, T ref is the reference temperature.

Impulsively started valve opening
In this section we consider the case where the valve is opened and the flow is impulsively started and compare the predicted velocity, temperature and Mach number for three different pressure inlet boundary conditions presenting three different cylinder pressure reserve at 75, 150 and 300 bar.
In Figs. 4-6, the maximum non-dimensionalised temperature, maximum non-dimensionalised mean velocity and maximum Mach number are presented, increasing the pressure in the cylinder from 75 to 150 and them to 300 bar has substantial effect on the fluid dynamics and fluid thermodynamics. As seen in Fig. 5, the fluid is surging due to the pressure difference and with the case of 300 bar the velocity of the gas reaches 2.7× the reference bulk velocity U b . The rapid increase of velocity is accompanied with an increase in temperature of 6× the reference temperature T ref , Fig. 4. As expected the simulation also exhibited that by decreasing the reserve pressure the pressure surge phenomena is decreased too. For both cases 300 and 150 bar the flow passed from subsonic to supersonic state crossing Ma = 1.0, this transition did not arise for the case of 75 bar, Fig. 6.
In Fig. 11(a)-(h) present across-sectional visualisation of the gas front as it advects through the valve's body. The temperature starts to rise as a result of the flow surging from the inlet toward the valve body, which can be seen in Fig. 4 as the first peak in the temperature profile and also in the velocity profile in Fig. 5. This first peak is a consequence of the fluid passing by a geometrical convergence regime resulting in a choked flow effect as seen in Fig. 11(a).
Later, the temperature profile exhibits a second even larger peak as the fluid starts expanding rapidly inside the valve's body as seen in Fig. 11(c)-(f). This peak is seen clearly in the flow profiles in Figs. 4-6. Finally a dip in the temperature occurs as the gas reaches the wall of the valve and advects toward the outlet, Fig. 11(g)-(h).
In Fig. 7 we present the instantaneous velocity magnitude profile at four locations inside the valve's body loc1 − 4 as illustrated in Fig. 3(a) and them with the reference solution provided from Rotarex S.A. For loc1 we notice a slight disagreement between our solution and Open-Foam solution, where the reference solution is predicting a higher peak in the velocity ≈ 2.5× where the hrSPH predicted only ≈ 1.5× the bulk flow. The reason can be explained by the fact that the outlet in the geometry used for the reference solution is slightly extended compared to the one we used in this study. Despite this disagreement we disclosed a good one for loc1 later in time. The loc2 profiles for both studies are in good agreement, this conclusion can be broaden for loc3 too, except that the hrSPH solution predicts a higher peak, and for loc4 the hrSPH predicts a larger dip in the velocity profile.

Time based valve opening
To replicate the opening of the valve, a time based pressure inlet boundary condition function is used, define as, where p inlet is the pressure inlet boundary condition at time t, p ref is the reference target pressure at the inlet, λ is the coefficient factor which determines the time needed to fully open the valve. We study three valve opening time dependent scenarios, scenario 1 the valve is fully opened at tU b /L v = 0.00342 and for scenario 2 at tU b /L v = 0.0342. The maximum temperature, velocity and Mach number profiles for the two time-dependent opening cases along with the impulsively  started case are represented in Figs. [8][9][10]. As expected the impulsively started case shows a higher peak in the temperature due to the sudden release of the reserved energy at the inlet, compared to the accumulated energy with time for the other cases, Fig. 8. Fig. 11(a)-(h) present a two dimensional cross section of the fluid front as it advects through the valve's body. The temperature starts to rise as a result of the flow surging from the inlet toward the valve body, which can be seen in Fig. 4 as the first peak in the temperature profile and also in the velocity profile in Fig. 5.

Flow structures
The three dimensional profile of the mean velocity is illustrated in Figs. 12(a-f), where the volume of the fluid is advecting through the body of valve. Fig. 12(d) shows the fluid at it largest expansion in the valve body with a dimensionless velocity of 2.5× the bulk one, this is also observed in Fig. 5. Fig. 13 shows the time averaged vorticity magnitude in the y-plane for the impulsively started scenario at p inlet = 300 bar. We observe that a region of increased vorticity exists in the centre of the outlet effectively showing a compact vortex core. A production of circulation is observed in the region of the converged geometry behind the inlet. The vortex core and the circulation are presented as 3-dimensional isosurface of the vorticity in Fig. 14.

Geometrical sensitivity analysis
As we analysed in the previous Section 7.1 the flow quantities (velocity and temperature) exhibit two peaks, the largest one is when the gas is at it's largest expansion inside the valve's body Figs. 4-6. We modified the design of the valve aiming at increasing the damping effect to the energy release. We are restricted with keeping the volume of the valves body, the inlet and outlet opening diameter and the length of the valve, hence the only region left to be modified is the geometrical convergence region.
In the modified valve design, Fig. 15(a) and (b), we extended this region toward the valve's wall which will have an effect in increasing the chocked flow effects, as seen in Fig. 17 represented with the first velocity peek.
The modified geometry design introduced a sharp corner, as a result the vorticity is increased by 14% compared to the initial valve design. We note here that, the vorticity increment can causes cavitation in the valve body on the long run, however as the modified geometry has a positive impact on the performance and safety of the oxygen respiratory valve, resulted in a lower peak temperature of only 2.4× the vorticity impact can be neglected for the time being, as cavitation is not the topic of this study.
Increasing the geometrical convergence region has a damping effect on how rapidly the gas is extending/surging later on, which can be clearly seen in the flow profiles in Figs. [16][17][18]. As seen, the effect on the gas dynamic is considerable since the modified design is now associated with an increase of ≈ 2.4× in the temperature compared to 6.0 × , and the velocity is increased by 1.6× compared to 2.7 × . From this modified valve design simulation we conclude that the developed flow characteristics will have a favourable effect on the valve functionality and safety.

Discussion and conclusions
Gas-dynamic pressure surges and adiabatic compression phenomena are generally hard to predict numerically. In this contribution, we simulate the dynamics of gas flow pressure surge and adiabatic compression in the fitted geometry of a generic of a respirator oxygen valve provided by Rotarex S.A. using the hybrid remeshed smoothed particle method.
We show that pressure reserve capacity has a considerable effect on the simulated flow fields (velocity, temperature), as the temperature could rise 6.0× the reference temperature.
Under the circumstances of valve malfunction due to manufacturing (e.g. particle contamination or grease contamination) issues, or user mistreatment of the valve in the filling facility since there is a possibility of connecting returned cylinders at full pressure to the filling manifold, or operating the valves in the wrong order or opening the valve quickly which can generate a shockwave. With aforementioned occurrences the developing flow characteristics can negatively impact the valve functioning and safety (fracture, ignition).
A modified valve design was developed during this study, by modifying the geometrical convergence region the new design exhibits a damping effect on the gas dynamics and flow characteristics, which has a favourable effect on the valve functionality and safety.
We believe that our results carry two important messages. First, the hrSPH method is able to predict the dynamics of the gas along with the pressure surge and adiabatic compression phenomena. Second, the reserved pressure has a dominating influence on the fluid flow in the valve, creating a modification of the valve mission-critical, which we illustrated by introducing a modified valve design. The modified valve geometry exhibited favourable effect on the developed flow characteristics.         1.8

Mach number
p inlet =300 bar, initial valve design p inlet =300 bar, modified valve design Fig. 18. The maximum Mach number profiles simulated using both the provided valve geometry and the new valve design, for the impulsively started scenario p inlet = 300 bar.