Pressure Waves due to Rapid Evaporation of Water Droplet in Liquid Lead Coolant

1 Ishlinsky Institute for Problems inMechanics of the Russian Academy of Sciences, Ave. Vernadskogo 101, Bldg 1,Moscow 119526, Russia 2Bauman Moscow State Technical University, Baumanskaya 2-ya 5, Moscow 105005, Russia 3National Research University “Moscow Power Engineering Institute”, Krasnokazarmennaya 14, Moscow 111250, Russia 4Joint Stock Company “Electrogorsk Research and Engineering Center on NPP Safety”, Saint Constantine Str. 6, Electrogorsk, Moscow Region 142530, Russia


Introduction
Lead-cooled fast reactors (LFRs) are considered a prospective avenue of nuclear industry development towards future systems possessing higher efficiency and operational and safety features in comparison with today's conventional light water reactors (LWRs) [1][2][3][4].Lead and lead alloys have several advantages over other liquid metal coolants (e.g., sodium) due to their chemically inert nature.Since no energetic reactions and hydrogen generation occur in molten lead coolant in the case of steam generator tube rupture (SGTR), reactor system design can be simplified significantly because the steam generator can be placed directly in the primary heat circuit, without the need of any additional heat transport systems.
The presence of high-pressure water in steam generator tubes raises concerns of possible consequences of tube rupture; therefore, extensive research has been focused on possible consequences of an SGTR accident.It is recognized that numerous physical phenomena accompany SGTR: shock wave propagation in the liquid lead; water flashing and expansion; water displacement from the tube to the intertube space filled with liquid lead; depressurization wave propagation in the broken tube, lead sloshing in the primary space; lead-water thermal interaction and possible steam explosions; steam bubble entrapment and transport towards the core [5].A number of dedicated experiments on waterlead interaction were carried out, most notably LIFUS experimental series [6,7]; these were supported by numerical simulations performed by SIMMER III thermohydraulic system code [8,9].It was found experimentally that water flashing produced significant pressure peaks in the liquid coolant (about 0.7-1.1 MPa reported in [7]); however, in SIMMER III calculations of these tests the initial pressure peak was underestimated 2-3 times (0.3-0.35 MPa).It should be noted that SIMMER III is a Eulerian multifluid code in which interactions at the scales of a single water droplet are not simulated directly but are taken into account via relevant closures.Such predictions can potentially be improved by more detailed (time and space-resolved) analysis of the phenomena involved in SGTR accidents.Similar experimental studies of SGTR accident in BREST reactor were carried out in Russia [10], with the main focus on the possibility of chain tube ruptures.The experimental results have shown that the chain ruptures of the neighboring tubes were unlikely, even if the tubes were initially defective.
The purpose of the current work is to consider an elementary process of SGTR accident, namely, rapid boil-up of a high-pressure water droplet immersed in low-pressure liquid lead.The aim is to obtain the features of two-phase mixture expansion, including the flow field and pressure waves generated by this expansion in the heavy liquid filling the surrounding space.First, a detailed model is presented requiring numerical integration of conservation laws in the water-vapor mixture and single-phase liquid lead, after which more simplified approaches are discussed relying on the solution of ordinary differential equations only.
Note that the processes involved are physically similar to those in the well-known boiling liquid expanding vapor explosions (BLEVEs) regarding one of the major hazards of accidents with high-pressure vessels [11].The difference, though, is that in the SGTR accident the ambient space is filled with a heavy liquid possessing much higher inertia and much lower compressibility than air.

Problem Formulation
2.1.Scheme of Multiphase Flow.Consider a spherical droplet of saturated water at a high pressure  ,0 and initial temperature equal to the saturation temperature,  ,0 =   ( ,0 ).At the initial instant, the droplet is immersed in a lowpressure space filled with a heavy liquid coolant (molted lead) at a lower pressure  ,0 (water is superheated with respect to the ambient conditions because  ,0 =   ( ,0 ) >   ( ,0 )).As the pressure relief wave is propagating through the high-pressure towards the droplet center, water boils up and expands due to rapid increase in the specific volume of the two-phase mixture, generating pressure waves in the ambient lead.In this work, we assume that the flow field is spherically symmetric and, therefore, consider the problem in the spherical coordinates with the origin at the droplet center.
The structure of the flow caused by expanding boiling water is sketched in Figure 1.Two main zones separated by a sharp interface are considered.Zone 1 corresponds to water in its single-phase (initial) state and in the form of twophase liquid-vapor mixture after the passage of boiling front converging towards the droplet center.Zone 2 corresponds to liquid coolant (lead).The expanding water works as a piston, causing propagation of a spherical pressure wave in the liquid lead.

Governing Equations.
In Zone 1, thermodynamic equilibrium between the liquid water and vapor is assumed.This assumption, discussed in detail in [12], implies that the interphase mass, momentum, and energy exchange processes are fast enough, so that both phases remain on the saturation line following the variation of local pressure; the velocity slip between the phases is also not taken into account (both phases have the same pressure and velocity).A further assumption is that the flow of the two-phase mixture is isentropic (i.e., at any point the sum of entropies of liquid water and vapor remains equal to the entropy in the initial all-liquid state corresponding to saturated conditions at the initial water pressure).
Under these assumptions, instead of considering the conservation laws for liquid and vapor phases coupled by the interphase exchange source terms, we consider the mixture of both phases characterized by its density   and velocity   for which the conservation equations are derived by summing up the corresponding conservation laws of both phases.In this derivation, the interphase exchange source terms (interphase drag, evaporation rate, etc.) entering the phase equations with opposite sign are cancelled out.As a result, the flow of boiling water in Zone 1 is governed by the Euler-type continuity and momentum conservation equations where  is the pressure,  is the time, and  is the radial coordinate.
The equation of state (EOS) for the water-vapor mixture is where V , and V ,V are the specific volumes of liquid and vapor phases on the saturation line depending on pressure .The mass fraction of vapor,  V , is determined as a function of local pressure  by the isentropic relationship where  ,0 is the specific entropy of saturated water at the initial water pressure  ,0 , while  , and  ,V are the specific entropies of saturated liquid water and vapor at the current pressure , respectively.Therefore, Note that  V and, therefore, all mixture properties depend parametrically on the initial pressure  ,0 ; otherwise all mixture properties are single-values functions of pressure , that is, the flow in the two-phase zone is barotropic.
For such a flow, the continuity and momentum equations ( 1) are decoupled from the energy equation.In fact, energy conservation was taken into account in the derivation of the equations of state ( 2)-( 4), rather than formulated as an independent differential equation.For example, such an important thermodynamic substance property as the heat of evaporation Δ V is contained implicitly in isentropic relationship (4): according to thermodynamics, the denominator (the difference of specific entropies of vapor and liquid on the saturation line) can be written as  ,V − , = Δ V /  , where   is the temperature of phase change (saturation temperature).
In Zone 2, liquid lead is considered as an inviscid compressible fluid; the flow is assumed to be barotropic (adiabatic) and governed by the continuity and momentum equations, similar to those in Zone 1 (see (1)): where ,   , and   are the pressure, density, and velocity of molten lead, respectively.For molten lead, the classic Tait equation of state (e.g., [13]) is used in the form where  * and  * are the pressure and density of lead at the reference conditions,  and  are EOS parameters specific to each substance.The Tait equation is widely used to describe compressible liquids in a wide range of pressure variation; for example, in underwater explosions, its application to molten lead is discussed below (Section 2.5.2).
To sum up, the mathematical model for high-pressure water droplet expansion in molten lead is formulated as Euler-type equations in both zones, but with different equations of state reflecting significant differences in the substance properties.In order to match the solutions in both zones, conditions on the moving interface between them must be formulated. 1 correspond to different substances (water and lead); they are divided by a moving sharp interface (contact discontinuity).The interface pressure,   , as well as its propagation velocity,   , is obtained as a solution to the Riemann problem.Given the instantaneous "left" and "right" states in water and lead, respectively, the following relationships determining the interface conditions can be obtained by the method of characteristics (see, e.g., [14,15]):

Interface Conditions. Zones 1 and 2 in Figure
where the top formula corresponds to water having pressure   and velocity   ; the bottom one corresponds to lead on the other side of the interface (with similar quantities denoted by subscript );  , and  , in ( 7) are the isentropic sound speeds in the respective substances.For the particular cases of EOS defined by ( 2)-( 4) and ( 6), the system of equations ( 7) for interface conditions takes the form (see [14]) where () =  −  ,0 +  (see (6)).The upper integral for water was evaluated numerically due to lack of analytical expressions for   () and  , ().Solution of equation set (8) gives the interface pressure,   , and its propagation velocity   .

Initial and Boundary Conditions.
The initial conditions correspond to a water droplet of radius  0 immersed in an infinite space filled with liquid lead; both fluids are assumed to be at rest: The interface is located initially at  =  0 , where pressure discontinuity exists, dividing saturated water and liquid lead.The boundary conditions are zero water velocity at the origin (due to spherical symmetry) and zero lead velocity "at infinity."The pressure in the origin is determined by the characteristic relationship (see the top formula in (7)), where the velocity on the left-hand side must be set to zero.At infinity, pressure tends to its ambient value in the liquid lead.

Substance Properties.
Consider in more detail the properties of each phase described by the respective equations of state: (2)-( 4) for water and ( 6) for liquid lead.

Saturated Water-Vapor Mixture.
In order to close the EOS for water, the pressure dependencies for the specific volumes, V , (), V ,V (), and entropies,  , (),  ,V (), of the liquid and vapor phases on the saturation line are required.These properties were taken from the NIST data [16] as appropriate tables within the pressure range 1-180 atm.Equally spaced pressure intervals of 0.1 atm were used, with linear interpolation of all properties between the tabular points.The same procedure was followed in [11] where explosions of pressure-liquefied gases in air were modeled.A possible alternative approach would be to implement any suitable correlations readily available for water (e.g., the well-known IAPWS formulations [17]).Consider now the specific features of a thermodynamic equilibrium two-phase mixture where both phases reside on the saturation curve, with the proportion determined by the isentropic condition.In Figures 2(a)-2(c), pressure dependencies of the mass fraction of vapor  V (), mixture density   (), and isentropic speed of sound  , () = (/  ) 1/2 =const are plotted for four initial pressures  ,0 in the range 10-20 MPa; see curves (1)-( 4) (note that since the initial state is assumed to be all-liquid, two-phase mixture exists only at pressures  <  ,0 ).Also, in Figure 2(c) the sound speed in the single-phase water and vapor is plotted by the dashed lines.One can see that density variation with pressure is quite considerable and nonlinear, the mixture possesses much higher compressibility than its constituent phases (mixture density variation with pressure mainly occurs due to phase transition, rather than compression of phases), and, therefore, the speed of sound in the mixture is much lower than in each phase.

Molten Lead.
Unlike the saturated water-vapor mixture (Section 2.5.1),thermodynamic parameters of molten lead for the conditions of interest in this work are far from its saturation line and, therefore, a single-phase formulation is appropriate.In this work, we use Tait EOS (6), which requires parameters  and  to be defined.No explicit data on these values was found for molten lead in the literature; therefore, the following approach was taken.The nominal density and speed of sound in molten lead at the nominal pressure  * = 0.1 MPa and initial temperature  0 = 800 K were evaluated from the experimental correlations [18,19]: The adiabatic speed of sound at the nominal conditions can also be obtained from Tait EOS (6): By equating the speed of sound ( 12) to ( 11) at the nominal conditions (  =  * ), we obtain that  =  *  2 , * =   = 3.212⋅ 10 4 MPa, where   is the adiabatic elastic modulus.Only the product of two constants can thus be derived from the density and speed of sound; therefore one of them must be set independently.Generally, this would require more extensive compressibility data than is available.However, in the current problem the molten lead pressure can only be as high as 18 MPa (the initial water pressure).Therefore, by recasting the Tait equation in the form 1/ (13) one can see that for typical values of  ≃ 10 (for water  ≃ 7) the second term in the brackets is actually very small, ( −  * )/ ≤ 5.6 ⋅ 10 −3 .Therefore, (13) can be presented in the linearized form   / * ≃ 1 + ( −  * )/() = 1 + ( −  * )/  independent of particular .This confirms that for the pressure range of interest in the current problem sensitivity to the parameters of the Tait equation is very low.
In the simulations we chose  = 12 (accounting for low compressibility of molten lead) and  = 2677.5MPa.

Numerical Method
The coupled system of equations ( 1)-( 6) was solved numerically.The one-dimensional computational domain of radius  max = 0.5 m, and a uniform mesh contained 1000 cells.In each subdomain, corresponding to Zones 1 and 2, relevant governing equations were discretized by an explicit highorder low-dissipation numerical scheme HR-SLAU2 belonging to the well-known AUSM family of numerical schemes [20].
In order to describe the sharp interface between Zones 1 and 2 (contact discontinuity), a conservative interface tracking method [14,15] developed for compressible multifluids was implemented.The interface was described by the levelset method [21]: a signed distance function was introduced, Another approach to taking into account the interaction on the contact discontinuity is the Ghost Fluid Method [22] which was also implemented and applied to test simulations.Both approaches gave very close results in terms of pressure profiles obtained.Note that in the current formulation (thermodynamically equilibrium water-vapor mixture) numerical procedure is free from the difficulties in approximation of stiff source terms describing the interphase processes (drag, evaporation, condensation, etc.) typically encountered in the multifluid models where conservation equations for each phase are solved separately.

Parameters.
As the baseline scenario, the parameters listed in Table 1 were selected.The initial pressures in water and coolant correspond to the operating conditions of Russian lead-cooled fast reactor BREST-OD-300 [3,4].The droplet size in the baseline scenario corresponds to the upper estimate of the volume of liquid water that can be discharged into liquid lead (10 −5 m 3 ) given in [5].It is comparable with the inner diameter of the steam generator pipes (13 mm) in this reactor.Keeping in mind significant uncertainties involved in such estimates, we also consider the case where a larger droplet is formed, of the radius 25 mm, which can be a result of a longer release time.Note that detailed study on water discharge process from a ruptured pipe would require a multidimensional model, which is beyond the scope of current study but will be performed in the future.

Pressure Waves Generated by an Expanding Water Droplet.
Consider the pressure and velocity fields developing during the first few milliseconds after droplet boiling and expansion begin.This stage is featured by quite fast processes of wave propagation from the high-pressure droplet into lowpressure lead, with the contact discontinuity on the waterlead boundary acting as a spherical piston.In Figures 3(a) and 3(b), radial pressure and velocity profiles are shown during the first 0.3 ms into the process.One can see that an elastic compression wave (generated by the initial pressure discontinuity) propagates through the liquid coolant at a high velocity, comparable with the speed of sound in liquid lead.This fast pressure wave has quite a sharp front; behind this front the velocity and pressure in the lead almost coincide at all times presented in Figure 3.This is further illustrated in Figure 4 where pressure time histories are shown at different distances from the droplet center: after the arrival of elastic compression wave, the pressure increases sharply and then remains almost constant.
Figure 5 shows in more detail the pressure profiles in the water droplet.It can be seen that pressure relief wave propagates into the saturated liquid converging to the droplet center at a time of about 0.1 ms, which means that by this time all water is boiling.As the boiling mixture expands into the surrounding coolant, the pressure in the droplet is falling gradually, with no significant radial gradients.In Figure 6, radial velocity profiles in are plotted near the expanding droplet.In the two-phase water mixture, the radial distribution of velocity is nearly linear, which means that expansion occurs more or less uniformly throughout the droplet.In the surrounding lead, the velocity profile is hyperbolic, as will be shown later.
A distinct feature of the processes considered, also mentioned in [5], is that, despite the high-pressure difference between the water droplet and surrounding liquid lead, the expansion velocity is rather low, reaching just few m/s.This is attributed to high inertia of the heavy liquid metal, in sharp contrast with explosive boil-up of superheated liquids in air (BLEVE phenomenon) where the interface expands at velocities comparable, or even exceeding the speed of sound in the air (hundreds of meters per second).This rather slow expansion is the main reason for the absence of significant pressure gradients in the two-phase water mixture, despite the low sound speed in it.This important finding can justify the application of a more simple approach based on the volume-averaged parameters for the water droplet.Simulations performed for a larger water droplet with the initial radius of  = 25 mm gave qualitatively similar results.Expectedly, amplitude of the pressure waves was higher due to larger energy yield for larger mass of water.In more detail, the pressure profiles will be compared in the following section where predictions will also be compared with a simpler model for water droplet expansion.

Comparison with a Simplified Model
Analysis of the results obtained above shows that, due to high inertia of the heavy liquid surrounding the water drop, its expansion proceeds much slower than propagation of the elastic compression waves through the molten lead and of pressure relief wave propagating through the twophase mixture in the boiling droplet.Water expansion can be considered as a relatively slow process with essentially subsonic characteristic velocities, superimposed by fast wave propagation at sonic velocities.Therefore, a simplified model can be derived for the "slow" expansion under assumptions that (i) in the two-phase mixture all parameters are distributed uniformly over the droplet volume, so that only the volume-averaged parameters of water-vapor mixture can be considered and (ii) molten lead behaves as an incompressible liquid.Such a model was recently developed in [23]; here we outline its main equations and focus on the comparison of the results with those obtained in Section 4 on the basis of complete problem formulation.
Under the above assumptions, the droplet is characterized by its radius  and mixture density   (  ) which depends on the current water pressure   according to the equation of state (2)-(4) (note that, unlike Sections 3 and 4, both quantities here are scalars rather than functions of radius).Water mass conservation implies that    3 = const, which gives a nonlinear relationship between the droplet radius and its pressure: Coolant is described as an incompressible liquid of density   = const, implying that in the ambient space ( ≥ ) the velocity   (, ) and pressure   (, ) satisfy the continuity and momentum equations By matching the pressure (  = ) and velocity (  | = = /) on the water-lead boundary and satisfying the boundary conditions "at infinity" (  = 0,   =  ,0 as  → ∞), we obtain a second-order ordinary differential equation of Rayleigh-Plesset type, widely applied for bubble oscillation and cavitation modeling [12]: Due to strong nonlinearity of   (  ) function (see Figures 2(b) and 2(c)), the system of coupled equations ( 14) and ( 16) must be solved numerically to provide the water droplet pressure   (), radius (), and expansion velocity () = /.After that, the velocity and pressure distributions in molten lead are obtained as In Figures 7(a) and 7(b), the pressure profiles obtained for two initial droplet sizes from the full model (see Sections 3 and 4) are compared with the corresponding solutions ( 17) and ( 18) at three time instants.The assumption of incompressible liquid adopted in the simplified model means that as soon as the interface starts to move, nonzero velocity   and nonconstant pressure   appear instantaneously in the whole space.The full model predicts finite-velocity propagation of compression wave; however, behind the wave front the pressure distribution practically coincides with that given by (18).Also note that during the short time interval for which results are plotted in Figure 7 the pressure profiles (18) are very close to each other, which means that the characteristic time scales of the processes described by the simplified model (radial outflow of coolant with no significant compressibility effects) are much longer than the characteristic times of compression wave propagation.
Thus, the full model (Section 3) must be applied in order to evaluate the effects of fast compression waves travelling in liquid lead at a speed of the order of speed of sound (about 1700 m/s).In such waves, sharp increase in flow pressure and velocity occurs; however, the duration of this transient impact is very short.Model ( 17), ( 18) is much simpler in implementation than the full model, and it can be used to study longer-term effects of the liquid lead flow caused by water droplet expansion.An example of such simulation is given in Figures 8(a) and 8(b) where time histories of the water droplet pressure and radius are shown for two initial radii, 13 and 25 mm (all other parameters were the same as those in Table 1).

Conclusions
The studies of single-droplet water-lead interaction show that effects on the surrounding structures (say, nearby tubes of steam generator) of high-pressure water expansion in molten lead can be of two types, differing in the characteristic time scales.The first type is a rapid (shock) impact caused by propagation of an elastic compression wave travelling at a speed close to the speed of sound in liquid coolant.Say, for a water droplet of 13 mm initial radius, the characteristic arrival time to a point located at 10 cm distance from the droplet is of the order of 0.05 ms, and pressure jump can be as high as 3.5 MPa.
Another type of interaction is the effect of slowly changing radial flow past the structures.This stage of interaction is featured by much longer characteristic times, of the order of water droplet oscillation time which can be as long as dozens of milliseconds.At this stage, liquid lead behaves as an incompressible fluid, and the force is mainly attributed to viscous drag.Evaluation and quantitative comparison of the effects of pressure waves and coolant flows on steam generator pipes will be the subject of future research.
While the 1D spherically symmetric statement of the problem used in the current work is the geometrically simplest, the proposed model takes into account the basic physical and fluid dynamics processes.In the future research, the model will be extended to more complicated 2D and 3D models in order to study various modes of water discharge from broken pipes into molten lead pool.

Figure 1 :
Figure 1: Sketch of the main zones for expansion of superheated water into liquid coolant.

Figure 5 :
Figure 5: Pressure profile near the expanding water droplet.

Figure 6 :
Figure 6: Radial velocity profiles at different times.

Figure 7 :
Figure 7: Comparison of pressure profiles predicted by the full (curves with times indicated) and simplified (curves denoted by symbols) models: (a) droplet radius  = 13 mm; (b)  = 25 mm.

Figure 8 :
Figure 8: Time histories of water droplet pressure, radius, and expansion velocity upon oscillations in liquid lead for different initial radii: (a)  = 13 mm; (b)  = 25 mm.

Table 1 :
Parameters of numerical simulation of single-drop waterlead interaction.withthe zero value reached on the interface.This function is transported on the mesh with the interface propagation velocity obtained from the solution of the Riemann problem (see Section 2.3), and approximations to numerical fluxes are modified appropriately in the cells cut by the interface.