An Experimental and Numerical Analysis on the Dynamical Behavior of a Safety Valve in the Case of Two-phase Non-flashing Flow

Pressure Relief Valves (PRVs) are key elements of any hydraulic system in the process industry, especially in chemical plants or hydraulic power transmission systems. Their task is to maintain the system pressure beneath a prescribed maximum pressure and vent the excessive fluid in an emergency scenario. This paper addresses the static and dynamic behavior of a Direct Spring-Operated PRV of conical shape in the presence of two-phase non-flashing flow, that is, water-air mixture. First, experimental results on the force and discharge characteristics of such a valve in a wide range of the air-to-water mass fraction are presented. Our test facility includes a custom-designed PRV with 42.5 mm inlet pipe diameter, an inlet pressure up to 6.6 bar(g) and a maximum lift of 10 mm. Additionally, the empirical results on the static characteristics, notably fluid force on the valve disc and discharge coefficients are reported as a function of the liquid mass fraction and valve lift. In the second part of the paper, we present the development of a Matlab-based simulation tool that is capable of predicting the dynamics and stability of such a valve in the case of two-phase, non-flashing, frozenmixture flow. Moreover, the effect of system parameters, such as spring stiffness and reservoir capacity are recorded. Finally, we also present results on the stability of the opening and closing the multi-phase flow influence on the stability of the blowdown process.


Introduction
Pressure Relief Valves are fundamental parts of any hydraulic system, which vent either single-phase flow (i.e. water, air), or, in more complicated cases, multi-phase flow (e.g. wet steam, water and air), as stated in [1][2][3][4][5][6][7]. The aim of such valves is to protect the system from excessive pressure by venting the unnecessary amount of flow if needed. However, sizing of these valves -especially in the case of multi-phase applications -might be challenging both from the static and dynamic point of view [8][9][10][11][12][13][14][15].
A PRV is an example of a single DoF oscillator which is prone to self-excited oscillators [16][17][18] and, if coupled with the fluid dynamics inside the piping system, might result in a surprisingly rich dynamic behavior, see [10] for more details. A significant research effort was devoted to predicting the dynamic behavior of safety valves in case of single-phase flow by applying CFD models or running empirical investigations. Despite the advantages of CFD tools to predict general flow behavior (see, e.g. [19]), the model setup and the computation requires a significant effort, not to mention parametric runs. The price of high-fidelity results of CFD comes at a price of scenario-based simulation and a lack of general understanding of the tendencies due to parameter variation. Hence, empirical and low-order models are needed for such a parametric analysis.
When trying to cope with the dynamical behavior of such valves, one needs to capture the fluid force and the amount of discharge. While the latter one is (relatively) straightforward to describe (as in [6,7,20]), different approaches were employed to model the fluid forces on the valve disc; the concept of "effective area" was introduced by Hős et al. [11] while other researchers use the flow deflection angles, as by Darby [12]. In both cases, single-phase cases were considered. This paper discusses the problem of predicting the static and the dynamic behavior of a PRV in the case of twophase, two-component flow with constant mass fractions (i.e. frozen mixture), for PRV design with the conical valve body. Even though we assume no interphase mass transfer, our model can also be applied for bubbly flows or the airflow contains a humidity [21][22][23][24] where the mass fraction remains constant. We assume that the force and discharge characteristics measured under steady state circumstances will remain valid (or give a reasonably good prediction) under dynamic scenarios as well [25]. We are particularly interested in the effect of the mass fraction on the stability of valve behavior.
Through the work, we employ DIER's ω technique to compute the theoretical mass flow and sonic velocity for a mixture composed of two components. This model is one of the HEMs models, which mainly due to the work of Leung [26]. The omega widely used to characterize a wide range of fluid conditions either the flashing flow (e.g. saturated liquid, sub-cooled liquid, and vapor-liquid mixture) or non-flashing (e.g. frozen two-phase), see [27][28][29] for further details. Omega method assumes both phases are in mechanical and thermodynamic equilibrium and neglects the velocity difference between the phases (no-slip assumptions). Moreover, it can be applied to both choked flow and non-choked flow. The only dependence of omega parameter lies on the fluid inlet conditions (i.e. stagnation conditions). Besides the advantages mentioned above, omega technique has some limitations such as losing the accuracy of the gained values at the critical point [30]. However, omega full details regarding the governed equations and the computations of the mass flux in the case of the frozen mixture can be found in [31]. This paper is organized as follows. Section 2 starts with the literature overview. Then, the SRVs description in Section 3 is introduced, while Section 4 presents the details of our test rig. In Section 5, the experimental (static) results are presented whereas Section 6 reports on the dynamical analysis. Finally, Section 7 summarized the investigation.

Literature overview
The protection and the reliable operation of any plant in process industries are the main goals of all manufacturers, and accordingly, the PRV is a cornerstone in any such system to achieve these purposes. Hence, the accuracy demand for Pressure Relief Valve sizing forms a necessity which in turn requires more experimental work, especially in case of multi-phase flows to find convenient and reliable approaches which help and enhance the aim of predicting the PRVs dynamics. Furthermore, the correct sizing ensures avoidance of unfavourable cases (oversizing or undersizing) which may occur during the blowdown of hydraulic systems and later may lead to shutdown of the system or even to the collapse of the facility (pipeline and tank). In contrast to the single-phase flow, most of two-phase flow researches are based on CFD models with few empirical attempts and in what follows we summarize some of the previous works to the author's knowledge with emphasizing the empirical works of main water to air mixture flow.
Dempster and Alshaikh [25] studied the two-phase (water/air) flow characteristics in the safety relief valve. Their experimental results indicate that the flow performance and disc forces are influenced by liquid mass fraction injection with a significant effect on the mass flow rate of the mixture. However, and in contrast to a previous limited data found in the literature, both authors found the mass fraction of liquid (water) has less impact on the disc forces.
Kourakos et al. [32] discussed the flow force in PRV for compressible, incompressible, and two-phase flow cases with focusing on valve opening positions. They compared the 2D model (CFD code) with experimental work, and they found a great agreement between the incompressible and the compressible flow results. However, the empirical results for air to water mixture under conditions of volumetric quality of air α = 20 % and set pressure (3 bar) at fully open position approach the compressible flow behavior. In contrast, the results between the proposed CFD model and the empirical outcomes show less matching as a result of measurement uncertainties.
Another experimental study for the critical flow of twophase flow (water/air) is examined in [33]. Dempster and Elmayyah [33] tested in this work an industrial refrigeration PRV over the pressure range of 6-15 bar(g), with water mass quality (mass fraction) (0 to 0.7), airflow (0.01-0.05 kg/s), and the discharge was near to the atmospheric conditions for different valve lift positions. They compared experimental results against a proposed CFD model. Moreover, the result of the CFD model also was examined against DIER's ω method and HNE-DS models adopted by the ISO standard [41]. The comparison shows that the mixture model (CFD model) flow rates have a good agreement with the experimental results at high valve opening positions for a low water mass fraction (0.11) with acceptable deviation (0.5 %). However, the accuracy of predicting the flow rate is less as liquid mass fraction increases with deviation (9 %) at water quality (0.55). On the other hand, for lower valve lifts the deviation tends to increase and becomes larger (around Dempster and Elmayyah [34] also discussed the prediction of discharge flow in non-flashing (water and air) of safety relief valve. The authors compared the critical mass flux, which is computed by omega method and CFD mixture model through SRV against an experimental work.
An industrial refrigeration SRV was tested over the pressure range of 6-15 bar(g), with air mass quality (mass fraction) (0.1 to 1). The discharge occurs to the atmospheric conditions for fully open positions. Both of the CFD based two-phase mixture flow and the HEM omega predict the same results within the mass fraction interval. Furthermore, CFD model and omega give a satisfactory prediction of the critical mass flow rate of the gas for liquid quality up to 0.4. However, for larger liquid mass fraction up to 0.9, error of 19 % has been recorded.
Arnulfo et al. [35] compared several experimental data from literature of two-phase flow and subcooled liquid of safety relief valve with different Nozzle models such as (HEM, ω-method, SEMs, Henry-Fauske's Homogeneous Non-Equilibrium, ERM, HFM, HDI, HNDI). They used the discharge coefficient from the models developed by Lenzing, Darby, and Leung. The comparison shows that for flashing flow, the model prediction by multiplying the Lenzing discharge coefficient with the theoretical mass flux estimated by HNDI (homogeneous non-equilibrium direct integration) resulting in good predictions. Similarly, applying the Leung/Darby formula regards the discharge coefficient, and considering the theoretical mass flow rate evaluated by omega method leads to the best prediction of non-flashing flow (water and air).
Burhani and Hős [36] addressed the effect of the frozen flow on the PRV dynamic behavior by employing DIERS' omega to predict the theoretical mass flow rate w.r.t different gas mass fraction ratios. The momentum force and the opening time predictions in term of effective area of three different geometries (two reaction-like valves with two different jet angles, and cone valve geometry) were discussed. They noticed an interesting opening time changes against very ridiculous gas mass fraction varying between few tenths of seconds up to few seconds. Additionally, the authors provided an analytical estimation of the valve opening time for both the existence and the absence of momentum force, which showed a great matching with simulated results with a maximum error of approx. 25 %. Fig. 1 depicts the cross-section of a typical DSO-PRV (Direct Spring-Operated Pressure Relief Valve) which comprises mainly of the body, Bonnet, and the movable components. The body of the valve holds all other parts, and the Bonnet is connected with the valve body by welding or threading to create a valve enclosure. The movable parts include the seat disc (to ensure leakage-free closure) and a disc holder loaded by an adjustable pre-compressed spring. Furthermore, spindle, spring, and a blowdown (closing/reseat pressure) adjustment ring to adjust the desired set pressure. If the pressure reaches the set point of pressure (which is adjusted by the pre-compressed of the spring), the valve opens and vents the undesirable flow rate of working medium towards the discharge side until pressure inside the hydraulic system drops to acceptable/safe levels (usually lower than the set pressure), then, the valve tends to close and to shut. The difference between the opening and reseating pressure is such so-called the blowdown, which is essential to predict high-frequency opening-closing cycles.

Equation of motion
The motion equation governs the PRV performance is a classic example of a single DOF oscillator; which is given by where m is the mass of the moving parts (disc, disc holder plus one-third of the spring mass), x x x , ,   are the valve lift, velocity, and acceleration respectively. k stands for the viscous damping coefficient, while s is the spring stiffness, and finally x 0 represents the pre-tension of the spring.
The F( x , p 0 , p b ) term refers to the total fluid force which acts on the valve disc, which combines three effects: 1. the force due to pressure distribution, 2. the force due to the momentum change of fluid (deflection) and 3. the force due to the wall shear distribution on the valve disc where p 0 is the upstream pressure and p b refers to the backpressure.
For simplicity, in this study, we assume that the dominant backpressure is the ambient pressure since the discharge of the fluid in our model (the designed geometry in Section 4) occurs towards the atmosphere (see [36]).
When the valve starts to open, the reaction force of the fluid on the valve disc can be defined by applying the momentum theory [37] with some assumptions which include the considering of a fixed and non-deformable control volume (Fig. 2) with neglecting the unsteady term (deformable volume) since it is has a small influence and these assumptions are supported by previous literature in [38,39], and hence we can write A v stands for the valve seat area, ρ in , ρ out , v in and v out are the fluid densities and velocities at the inlet and the outlet of the control volume respectively and θ is the deflection angle (jet angle).
To give more generality for the findings, and to show the effect of the disc forces, results are provided in term of dimensionless approach [40]; we define the effective area (force coefficient) and hence Eq. (2) can be written where the effective area is The effective force is straightforward to measure as it requires only the fluid force and the pressures at the upstream and discharges side. Generally, the effective area curve depends on the valve lift primarily, whereas, for compressible fluids, the density change and the presence of choked flow condition introduce further (secondary) pressure dependence.

Test rig description and procedure
A hydraulic facility was built up in Budapest University of Technology and Economics at the laboratory of the Department of Hydrodynamic Systems to test a cone valve performance in the case of two-phase (non-flashing) flow. The measurements were performed at different valve displacement positions, and both the flow rates and the forces on the valve disc were measured for a wide range of waterto-air mass fraction (6.5 % to 30 %) under a steady-state condition with non-heat transfer assumption (short nozzle). The strategy of the current experimental work is to inject different water mass flow rates in the upstream of pressurized airflow to form a mixture of two-phase flow which later hits the valve disc. Fig. 3 (schematic of test rig) and Fig. 4 (real test rig) show the main parts of the test rig which essentially consists of the pipeline with an inner diameter D p = 42.5 mm and the length of 145 cm.
The pipe is also connected via an elbow to a large air reservoir with a capacity of V = 900 liter where the pressurized air comes from a compressor which has the properties of Q = 160 m 3 /hr, P = 22 KW and n = 980 rpm. Additionally, a Bourdon pressure gauge is attached to the tank to show the working pressure where the tank temperature (t = 25 °C) Along the main pipeline, a sonic nozzle with a throat diameter (d n = 10 mm) and two ports for absolute pressure transducers can be seen where the first pressure sensor is Furthermore, the water injection process can be achieved by utilizing one of two liquid spray nozzles (490.403CA-490.683CC LECHLER) with different capacities 1 and 5 liter/min respectively at 5 bar pressure difference (nominal pressure). The desired water mass fraction ratio, that is x m m m l l l = + ( )    air (subscript "l" stands for liquid), where x l + x air = 1, and   m m l , air , are the water and the air mass flow rates respectively. The two nozzles provide a uniform full cone axial-flow with spray angle 45 °C. The spray nozzle is positioned in the centerline of the pipe before the valve, while a plastic water hose feeds the nozzle. Also, a regulator and Bourdon pressure gauge are mounted on the hose to control the water flow rate and read out the water pressure magnitude as well. Likewise, the force can be measured by using force sensor KM300 made of Aluminum where the measuring range up to 12 kN. Last but not least, an MX840A HBM Data acquisition with eight channels is used to collect the measured data, with Smart Catman V3.5.1 version.
The measured data in this paper was taken at different valve lifts by using a movable mechanism (see Fig. 5) to set the desirable opening position, where those opening positions are corresponding to specific water mass fraction varies between pure air at x l = 0 and water to air mixture at x l = 30.8 %. Furthermore, the valve encounters airflow rates  m air kg/sec = − 0 11 0 12 . . for upstream pressure range p 0 = 5.5 − 6.6 bar(g). The temperature change of the airflow during the expansion via the sonic nozzle is neglected. The displacement range of the studying in this work is until we reach the critical valve lift at x = D p / 4, at which the flowthrough area of the valve becomes larger than the pipe area, hence the location of choking changes. Each measurement lasted for approximately ten seconds (between the valve opening and the valve closure) and we have cut out (manually) an internal interval of a few seconds with stabilities pressure levels for post-processing purposes.

.1 Effective area
The flow force (the first static characteristic) of five different liquid injection processes are shown in Fig. 6 under the definition of effective area term against variation in the disc lift x (the maximum valve lift is 10 mm) in term of relative valve opening which is  x x D p = 4 . The measured data is represented within the experimental error (error bar). Furthermore, the mentioned figure shows both the measured data points and the curve fits for constant mass fractions. We see that higher water mass fractions result in slightly higher forces, but the trends are the same A eff starts from a value near unity (the momentum force is almost zero, and the pressure force dominates) at small relative valve lifts  x < 0 1 . and drastically increases at  x = 0 2 . relative opening position before it decreases again until valve reaches  x = 0 8 . . The behavior of the curves are close to each other up to approx. 0.6 relative lift, especially for the first three liquid mass fractions. However, for the ratios x l = 17.2 %, 20.7 %, and 30.8 % the effective area shows higher values at the last third of the relative valve opening interval and even shows a slight increase. However, discovering the root cause of the behavior as mentioned earlier by using CFD tool would be exciting, but this is beyond the scope of the current work.

Discharge coefficient
The discharge coefficient is simply the ratio of the actual (measured) flow rate and the theoretical prediction assuming ideal flow. In this paper, the actual mass flow rate is evaluated by measuring the water flow by pressure measurement (by pre-calibrating the spray nozzle characteristic) and the air flow by the sonic nozzle, while the theoretical flow rate is computed by applying DIER's ω technique in case of non-flashing flow, that is (see [29] for details) where α 0 is the volume fraction of the gas phase (air) and can be given with the following expression (see [41]) where x g is the mass fraction of the gas, ρ g and ρ l are gas and liquid densities, respectively. Fig. 7 depicts our results within the experimental error bar. The discharge coefficient is decreasing with increasing the relative valve opening in the majority of the studied displacement interval, up to approx. 0.6. Beyond this relative lift, it seems to stabilize at approx. C d ≈ 0.5 for x l = 0 %, 6.5 %, 11.1 % and 17.2 %. On the other hand, the discharge coefficient values for the rest of the liquid mass fractions is approximately 0.6. This result is not surprising since the pressure range at the last third of the relative displacement range is close to the magnitude of ambient pressure (constant pressure).
We have to add that there are several weakness in our measurements. First of all, the displacement measurement is inaccurate for small valve openings and non-concentric valve displacement (relative to the pipe axis) might cause issues. Indeed, we have measured discharge coefficients slightly above 1, which were trimmed artificially. On the other hand, for large openings, the pressure drop across the valve is small, which causes again measurement inaccuracies. To highlight these issues, we have also added the estimated error bars to the data points.  Fig. 6 Effective area vs. relative valve lift Fig. 7 Discharge coefficient for liquid (water) to air mixture (ode45 solver), which is a 4th-order Runge-Kutta solver with built-in step size control. We have coupled the reservoir dynamics where V is the reservoir size,  m in is the inlet mass flow rate, C d is the discharge coefficient (see Fig. 7), A out (x) is the flow-through area and G( p 0 , p b ) and ( a s ) are the mass flux and sonic velocity of the mixture respectively, both of them estimated by DIER's ω technique. We have used the curve fits depicted in Fig. 6 and Fig. 7 to handle the effect of varying mass fraction. Table 1 shows the valve properties of the PRV model used in the simulation.
Typically, we have run two simulations: once neglecting the momentum force, that is, setting A eff = 1 in Eq. (3), these are the red time histories in both Fig. 8 and Fig. 9. A second run followed then with the measured A eff curves.

Effect of reservoir volume
As a first step, we have run several simulations with fixed spring stiffness but varying reservoir volume. Intuitively, the larger the reservoir volume is, the less important its dynamics will be as the pressure fluctuations inside the reservoir will decrease as its volume is increased. Some examples runs are depicted in Fig. 8 and Fig. 9. Fig. 8 and Fig. 9 depict the same venting scenario, but with different vessel sizes. The blue line in the bottom-most panels shows the (constant) inlet mass flow rate (into the reservoir), see Eq. (7). The red line in each figure depicts the simulation results with a constant effective area A eff = 1 assumption (see Eq. (1) and Eq. (3)), that is, the case when the force due to the momentum change on the valve disc is neglected. The black lines correspond to the case when the measured effective area curve was used, that is the curve-fit to the pink triangles of Fig. 6.
The middle panels depict the vessel pressure (the blue dash-dot line is the set pressure) while the top-most panels show the valve displacement.
The first difference we observe is the change in the valve lift (and vessel pressure) time history when using a small vessel (Fig. 8, V = 100 liter) and a large vessel (Fig. 9, V = 900 liter). It is clear that the larger the vessel is, the longer it takes the pressure to build up until it reaches the set pressure, but it is also interesting that in both cases, the contribution of the impulse forces (black vs. red lines) is fundamental. We observe a fixed valve opening in the case of no momentum forces (red curves, A eff = 1), while, in the case of realistic force estimation (black curves), we see a periodic opening-closing of the valve. Moreover, we also see that the opening pressure and reseat pressure (when the valve closes again), are different; in Fig. 8, we see that the reseat pressure is about 96.5 % of the opening (set) pressure. This well-known difference between the set and reseat pressure is called blowdown.   (See text for details.) Fig. 10 depicts the maximum valve lift ( Fig. 10 (a)) and the minimum and maximum pressure values ( Fig. 10 (b)), relative to the set pressure p set = 6 bar(a) as a function of the reservoir volume. Circles stand for the case without the momentum term while crosses denote the simulations performed with the measured effective area curves. The difference between the minimum and maximum pressure is essentially the blowdown.
For large enough vessels, neither the maximum lift nor the blowdown varies. However, if the vessel size is smaller than approx. 300 liters, we experienced instability even in the case constant effective area curve (see the increase in max(x) in the data points with circle). For larger vessels, we see a constant lift (over the time) in the case of the simulations without the momentum force and a periodic opening-closing in the case of the realistic effective area curve, which is referred to as cycling in the literature.

Effect of spring stiffness
As a next step, we have run several simulations with fixed reservoir volume but varying spring stiffness, adjusting the spring pre-compression x 0 in such a way that the set pressure remained constant. We again set off with the scenario depicted in Fig. 7, but now, while keeping the reservoir volume constant, we increased the spring stiffness up to a value of 20 times. The results are depicted in Fig. 11. The color coding is the same as in Fig. 6 and Fig. 7; the circles denote the cases without the momentum force (constant effective area) while the crosses denote the case with the measured effective area curves.
Notice that softer springs result in larger valve openings and instability. This behavior is well-known from the literature, see [8,9] for an example. It is also interesting to see that the gas case (black line) is somewhat more stable than the mixtures; the higher the water content is, the larger the unstable region and the amplitude of the oscillations become.
The lower panel highlights the blowdown behavior, which is the difference between the set pressure (due to the rescaling of p by p/p set , 1 in the figure) and the minimum pressure (see Fig. 8). The fact that in the case of no momentum force (A eff = 1, circles) there is no blowdown, while if the momentum force is considered, blowdown emerges (that is, the minimum pressure is lower than the set pressure) highlights the importance of the use of realistic fluid force model. It is also interesting to see that higher spring stiffness results in lower blowdown and, with high enough spring stiffness, blowdown is completely lost.

Conclusion
This study unveiled the effect of mass fraction of non-flashing multi-phase flow on the dynamics of Direct Spring-Operated Pressure Relief Valves. We have employed DIER's ω technique to predict the density, sonic velocity and mass flux of the flow.
Firstly, by means of experiments, we have shown that the mass fraction plays an important role in both the fluid force and the mass flow rate through a relief valve. Even though the tendencies are the same, one cannot omit the effect of impulse forces when trying to cope with dynamic scenarios such as opening and closing of the PRV. We have highlighted these effects by comparing dynamic simulations with and without the fluid impulse forces, e.g. Figs. 8-11. Secondly, we have shown that the reservoir size does not play a significant role, provided that exceeds a critical minimum value. On the contrary, the PRV spring stiffness plays a key role; soft springs result in self-excited oscillations while stiff springs eliminate the blowdown behavior of the valve.
Finally, by using DIER's ω technique, we have shown that is straightforward to include multi-phase effects in PRV simulations as long as the mixture can be handled as a (thermodynamically and mechanically) homogeneous phase. Our future plans include running simulations with more complicated multi-phase fluids, e.g. wet steam.