# Modeling of Spray System Operation under Hydrogen and Steam Emissions in NPP Containment during Severe Accident

Authors: Vadim E. Seleznev; Vladimir V. Aleshin; Sergey N. Pryalov
Journal: Science and Technology of Nuclear Installations
 (2011)

---

## Abstract

The paper describes one of the variants of mathematical models of a fluid dynamics process inside the containment, which occurs in the conditions of operation of spray systems in severe accidents at nuclear power plant. The source of emergency emissions in this case is the leak of the coolant or rupture at full cross-section of the main circulating pipeline in a reactor building. Leak or rupture characteristics define the localization and the temporal law of functioning of a source of emergency emission (or accrued operating) of warmed up hydrogen and steam in the containment. Operation of this source at the course of analyzed accident models should be described by the assignment of the relevant Dirichlet boundary conditions. Functioning of the passive autocatalytic recombiners of hydrogen is described in the form of the complex Newton boundary conditions.


---

## Body

## 1. Introduction

The paper describes an approach to fluid dynamics simulations of processes occurring inside the nuclear power plant (NPP) containment during operation of spray systems in severe accidents. This topic is relevant at the current stage of advancement in high-accuracy numerical simulations of NPP operation [1–3]. In the case of interest, accidental releases mostly occur as a result of a coolant leak or line rupture in the main coolant circuit of the reactor building. Leak or rupture characteristics determine the location and the time law of accidental release (or buildup) of heated hydrogen and steam into containment rooms. The development of the mentioned time law proper is beyond the scope of this paper. It is reasonable to simulate the behavior of a source of hydrogen and steam during the severe accident under consideration by specifying the relevant Dirichlet boundary conditions. Containment rooms are equipped with passive autocatalytic recombiners (PARs) of hydrogen and condensation heat exchangers of the reactor building’s passive heat removal system. PAR modeling techniques have been discussed in detail in [4–6]. In our case, PAR operation is described using the complex Newton boundary conditions, the time history of which is defined using the known submodeling principle [4–6]. Condensation heat exchangers are modeled by specifying the Newton boundary conditions, the algorithm of generation of which is basically similar to the definition of heat exchange conditions on inner containment walls accounting for their steel liner (see Section 5).The models proposed in the paper are intended for engineers and technical staff involved in the NPP design and maintenance at both construction and operation stages. Numerical analysis of these models can be performed without high-performance computers.One of the mandatory prerequisites for the mathematical modeling in our case is that we should correctly describe the operation of the spray system and heat transfer in steel-liner containment walls accounting for steam condensation on them. The spray system is installed under the dome of the central reactor building. According to its purpose, it belongs to confinement safety systems, and its basic functions include pressure and temperature control in the containment during severe accidents.Combustion risk of such clouds in this case is assessed using the known Shapiro diagram. In the first approximation, the aqueous solution of boric acid distributed by the spray system can be considered similar to water in its physical properties.Within our model we assume that the gas mixture in containment rooms has six components: hydrogen produced during severe accidents,(m=1); air as a mixture of its three basic components—oxygen, (m=2), nitrogen, (m=3), and atmospheric water vapor, (m=4); steam injected into the containment during severe accidents, (m=5); additional steam produced by hydrogen and oxygen absorption by PAR [4–6] (referred to as absorption product) (m=6).


## 2. Selection of Underlying Model Equations

As a basis for the development of the fluid dynamics model of processes inside the containment, we used an adapted set of Reynolds equations completed by the(k-ω)-turbulence model [7–9].The turbulence model was chosen by the authors based on the results of computational experiments with steam-hydrogen-air flows in containment rooms without considering the operation of the spray system. These experiments demonstrated considerable advantages of this model over the zero-order turbulence models and various modifications of the(k-ε)-turbulence model [10, 11]. Similar conclusions have been made in [12, 13]. Application of the zero-order turbulence models in most of computational experiments lead to clearly physically invalid numerical results for the space-time field of gas flow velocities in the simulated rooms.Thus, the basis for the fluid dynamics model of processes inside the containment with neglect of the spray system operation can be represented in the following way (see, e.g., [7–9, 14]):(1*)∂ρ∂t+∇⃗⋅(ρV⃗)=0,(2)∂(ρYH2)∂t+∇⃗⋅(ρYH2V⃗)=∇⃗⋅(μ+ρνTScH2∇⃗YH2)-ωH2,(3)∂(ρYO2)∂t+∇⃗⋅(ρYO2V⃗)=∇⃗⋅(μ+ρνTScO2∇⃗YO2)-ξO2massωH2,(4)∂(ρYN2)∂t+∇⃗⋅(ρYN2V⃗)=∇⃗⋅(μ+ρνTScN2∇⃗YN2),(5)∂(ρYatm_H2O)∂t+∇⃗⋅(ρYatm_H2OV⃗)=∇⃗⋅(μ+ρνTScatm_H2O∇⃗Yatm_H2O),(6*)∂(ρYH2O)∂t+∇⃗⋅(ρYH2OV⃗)=∇⃗⋅(μ+ρνTScH2O∇⃗YH2O),(7)YPAR_prod=1-YH2-YO2-YN2-Yatm_H2O-YH2O=1-∑m=1N-1Ym,(8*)∂(ρV⃗)∂t+∇⃗⋅(ρV⃗V⃗)=(ρ-ρatm)g⃗-∇⃗P+∇⃗⋅(μ+ρνTμτ)-23∇⃗(ρK),(9*)∂(ρH)∂t+∇⃗⋅(ρHV⃗)=∂P∂t+∂Q∂t+ρg⃗⋅V⃗+∇⃗⋅[(λ+λT)∇⃗T+μ+ρνTμτ⋅V⃗-23ρKV⃗]+∑m=1N∇⃗⋅[(Cp)mTμ+ρνTScm∇⃗Ym],(10)∂(ρK)∂t+∇⃗⋅(ρKV⃗)=∇⃗⋅[μ+ρνTPrK∇⃗K]+G-β*ρKω-νTPrTg⃗⋅∇⃗ρ,(11)∂(ρω)∂t+∇⃗⋅(ρωV⃗)=∇⃗⋅[μ+ρνTPrω∇⃗ω]+ωK[νTPrTαG-β**ρKω+νTPrT×{(α+1)max⁡(g⃗⋅∇⃗ρ;0)-g⃗⋅∇⃗ρ}νTPrT],(12)P=ρR0T∑m=1NYmMm+ρatmg(x3-x3,0),H=CpT+0.5V⃗⋅V⃗;α̃S=∑m=1Nfα̃(Ym,α̃m),α̃≡μ,λ,Cp,(13)μ=(T273.15)3/2273.15+CST+CSμ0,λT=CpρνTPrT;Scm=μρDm,(14)τij=μ[(∂ui∂xj+∂uj∂xi)-23δij∂uk∂xk],(15)G=ρνT{12(∂ui∂xj+∂uj∂xi)2-23(∂uk∂xk)2}-23ρK∂uk∂xk,νT=α1Kmax⁡{α1ω;S̃F2},S̃=2S̃ijS̃ij,S̃ij=12(∂ui∂xj+∂uj∂xi),F2=tanh[max2{2Kβ*ωy;500vy2ω}],
where D(⋯)/Dt≡(∂(⋯)/∂t)+V⃗·∇⃗(⋯) is the substantial derivative of the scalar function (writing a substantial derivative of a vector function means substantial differentiation of vector function components); t is time; V⃗ is velocity with components u1, u2, u3; ∇⃗ is nabla; ρ, ρm, ρatm are density of mixture, mth mixture component, and atmosphere at some point in space of containment rooms, respectively; YH2=Y1, YO2=Y2, YN2=Y3, Yatm_H2O=Y4, YH2O=Y5, YPAR_prod=Y6, Ym=ρm/ρ is the relative mass fraction of the mth component; μ is the dynamic viscosity; μT=ρνT is the turbulent viscosity; ν is the kinematic viscosity; Scm is the Schmidt number corresponding to the mth component; ωH2 is the mass rate of the generalized irreversible single-stage reaction of hydrogen absorption in PAR [4, 5]; ξO2mass is the stoichiometric mass coefficient of oxygen; N is the number of components in the gas mixture (in our case, N=6 (see above)); the subscript T stands for “turbulent,” and the subscript atm, for the “atmosphere in containment rooms”; g⃗ is gravity acceleration (g is the gravity acceleration modulus); P is piezometric pressure of mixture; τ is tensor of viscous stresses with components τij; K is turbulent kinetic energy; H=h+0.5V⃗·V⃗ is total specific (per unit mass) enthalpy (h=CpT is enthalpy for perfect gas, Cp is specific (per unit mass) heat capacity at constant pressure; T is mixture temperature); ∂Q/∂t is the specified rate of heat release from outer sources per unit volume (heat generation rate per unit volume); λ is the thermal conductivity; λT is the turbulent thermal conductivity; Pr is the Prandtl number (the subscript K indicates that the value of the Prandtl number is defined specifically for turbulent energy equation (10), and ω indicates that it is defined for turbulence dissipation (11));  PrK=2.0 [7–11]; β*=0.09 [7–9]; ω are turbulence frequencies; Prω=2 [7–9]; α=5/9 [7–9]; β**=0.075 [7–9]; R0 is the universal gas constant; Mm is molar mass of the mth component; x3 is the vertical coordinate in the Cartesian frame of reference directed away from the earth center   (here, x1, x2, x3 are radius vector coordinates of the point in the Cartesian frame of reference); x3,0 is the coordinate of the point x3 with a defined value of ρatm;  fα̃(⋯) are known semiempirical functions; CS is the Sutherland constant; μ0 is the dynamic viscosity under normal conditions; Dm is the coefficient of binary diffusion of the mth component in the rest of the mixture; δij is the Kronecker delta; α1=5/9 [7–9]; y is the distance to the nearest wall. The functions G, S̃, F2 in (15) are auxiliary. Note that the model uses Favre-averaged velocity components and thermal variables (H, h, T) [15] (i.e., mixture density is used as a weighting function) and classical Reynolds-averaged density and pressure [11]. Asterisks in the right upper corner in (1)–(15) and in the following indicate the equations that will be modified as the model will be constructed.It should be noted that in the model (1)–(15) steam was conventionally divided into 2 components, such as atmospheric steam and injected (during accident) steam. That was done for the implementation of comprehensive analysis of distribution parameters of injected (during accident) mixture in containment rooms.


## 3. Modeling of Spray System Operation

In order to account for heat exchange processes during production (evaporation) and condensation of steam by the spray injection, model (1)–(15) needs to be modified. The first step of such modification will include converting equations (1*), (6*), and (9*) to the following form: (1)∂ρ∂t+∇⃗⋅(ρV⃗)=ωH2O;(6)∂(ρYH2O)∂t+∇⃗⋅(ρYH2OV⃗)=∇⃗⋅(μ+ρνTScH2O∇⃗YH2O)+ωH2O;(9**)∂(ρH)∂t+∇⃗⋅(ρHV⃗)=∂P∂t+∂Q∂t+0.5Qphase_transition×[ωH2O-|ωH2O|]+(Cp)H2OTsat⁡ωH2O+ρg⃗⋅V⃗+∇⃗⋅[(λ+λT)∇⃗T+μ+ρνTμτ⋅V⃗-23ρKV⃗]+∑m=1N∇⃗⋅[(Cp)mTμ+ρνTScm∇⃗Ym],
where ωH2O≥0 is the mass rate of steam production as a result of evaporation of water droplets distributed into the reactor building by the spray system; ωH2O<0 is the mass rate of steam condensation on water droplets distributed into the reactor building by the spray system; Qphase_transition is specified latent mass heat of steam condensation into water droplets at ωH2O<0; Tsat is temperature of saturated steam. In fact, the quantity ωH2O is the mass of steam produced (or condensed) per unit time in unit volume of mixture.In our case, evaporation of each droplet is caused by the heat supplied to it from the surrounding gas mixture and spent on the phase change due to the given latent mass heat of water droplet evaporationQphase_transition at ωH2O≥0. There is no heat transfer from the droplet to the surrounding gas mixture. To include this, (9**) uses the summand 0.5Qphase_transition[ωH2O-|ωH2O|]. The summand in (9**) given by [(Cp)H2OTsatωH2O] describes enthalpy increment in the gas mixture due to the steam newly produced from droplets or enthalpy decrement with steam disappearance as a result of its condensation on droplets. Note that the heat release of condensation of superheated steam is somewhat higher than that of saturated steam. This difference, however, does not exceed 3 percent [16], so it can be reasonably ignored in applied simulations.We obtain a formula for the rateωH2O with droplet evaporation (ωH2O≥0). [17] established the equations that characterize the rate of steam production J from a spherical droplet: (16)J=J1=ρDsteam0.5ddropln⁡[1-Ysteam,∞1-Ysteam,sat⁡],(17)J=J2=λsteam0.5ddrop(Cp)H2Oln⁡[1+(Cp)H2O(T∞-Tsat⁡)Qphase_transition],
where ρ is gas mixture density; Dsteam is the coefficient of steam diffusion in the mixture; ddrop is droplet diameter; Ysteam,∞ and T∞ are the mass fraction and temperature of steam far from the surface, respectively (where effects produced by the evaporating steam can be ignored); Ysteam,sat is the mass fraction of saturated steam; λsteam is the thermal conductivity of the gas (steam) mixture. Note that (17) can be extended only to the processes of droplet evaporation, while (16) is applicable to both evaporation and condensation processes (see below). We assume conventionally that all droplets in the small finite volume are of the same diameter. The droplet surface area is Sdrop=πddrop2,  where π is Pythagorean number. Consequently, using formulas (16) and (17), we can estimate the mass of steam produced from a single droplet per unit time:(18)dmsteam_dropdt=-dmdropdt=πddrop2J,
where dmsteam_drop/dt is the rate of steam production from one droplet; mdrop is droplet mass; dmdrop/dt is the rate of droplet mass variation. Let n be a specific (per unit volume of matter) number of spherical water droplets. In this case, the mass of steam produced from all unit-volume droplets can be calculated using the formula (see (17)) (19)ωH2O=ndmsteam_dropdt=-ndmdropdt=nπddrop2J=2nπddropλsteam(Cp)H2Oln⁡[1+(Cp)H2O(T-Tsat⁡)Qphase_transition].
Here, as with (1)–(15), the temperature T is equal to the gas mixture temperature (far from the droplet surface).As it evaporates, the droplet decreases in its diameter. Let us estimate the rate of droplet diameter variation. The droplet mass is(20)mdrop=ρliquidVdrop=ρliquidπddrop36,
where ρliquid is density of the droplet liquid (water); Vdrop is droplet volume. Hence,ddrop=[6mdrop/(πρliquid)]1/3, and
(21)ddtddrop=13[6πρliquid]1/3mdrop-2/3⋅dmdropdt=-ωH2O3n[6πρliquid]1/3mdrop-2/3.
These formulas allow us to analyze the process of droplet size variation with evaporation. This process occurs, if the pressure of saturated steam near the droplet exceeds the steam pressure in the gas mixture. Proposed model does not take into account the effect of droplet collision and the interaction between water droplet and wall.However, at the initial stage of its flight, the droplet has a low temperature, at which the steam pressure can be lower than the steam pressure in the gas mixture. For this reason, during the initial period of the flight, steam can condense on the droplet surface with droplet warming up due to latent heat of steam condensation. Let us consider the process of steam condensation onto the droplet(ωH2O<0). Because of the small size of droplets distributed by the spray system, the temperature of each droplet is assumed to have its own constant value throughout the droplet volume, including its surface. As a result, the saturation temperature (just as the gas mixture temperature) near the droplet is equal to the droplet temperature: Tsat≈Tdrop. Let us derive an equation to describe the droplet temperature variation during steam condensation on the droplet surface. The liquid is assumed to be governed by the ideal caloric equation of state (EOS), hliquid=ρliquidCliquidTliquid, where hliquid is specific (per unit volume) enthalpy of the liquid; ρliquid=const is liquid density (the liquid is incompressible); Cliquid=const is specific (per unit mass) heat capacity of the liquid; Tliquid is liquid temperature. Considering this EOS, droplet enthalpy can be written as Hdrop=mdropCliquidTdrop=ρliquidVdropCliquidTdrop, where Hdrop and Tdrop are average enthalpy and temperature of the droplet, respectively. Hence, (22)ddtHdrop=Cliquid(mdropdTdropdt+Tdropdmdropdt)orddtHdrop=Cliquid(mdropdTdropdt-TdropωH2On).
On the other hand, condensation is accompanied by heat gain due to steam condensation and growth of droplet mass due to the condensate with the saturation temperature Tsat:(23)ddtHdrop=(Qphase_transition+CliquidTsat⁡)dmdropdt=-(Qphase_transition+CliquidTsat⁡)ωH2On.
By equating the right-hand sides of the second equation in (22) and (23), we have(24)dTdropdt=ωH2Onmdrop[(Tdrop-Tsat⁡)-Qphase_transitionCliquid]=-ωH2OQphase_transitionnmdropCliquid.Now, let us derive an equation to calculate the mass rate of steam condensationωH2O on water droplets (ωH2O<0). Formula (16) results from solving the component continuity equation for a steady-state gas mixture flow(25)∇⋅(ρwYsteam-ρDsteam∇Ysteam)=0
for the following boundary conditions: r=r0=0.5ddrop:Ysteam=Ysteam,sat,  r=∞:Ysteam=Ysteam,∞, where r is distance from the droplet center to an arbitrary point.During droplet evaporation,Ysteam,∞≤Ysteam,sat, which produces a flow of steam from the droplet surface to the gas mixture. If the situation is opposite, Ysteam,∞>Ysteam,sat (i.e., the steam mass fraction in the gas mixture is higher than on the droplet surface equal to Ysteam,sat), a reverse flow of steam—from the gas mixture to the droplet surface—should occur. The continuity equation (25) for the mixture component will remain valid. Consequently, in order to calculate the mass rate of steam condensation ωH2O on water droplets, one can use (16) (see also (19)). Hence, (26)ωH2O=nπddrop2J=2nπddropρDsteamln⁡[1-Ysteam,∞1-Ysteam,sat⁡].
If we first write the mass fraction of steam on the saturation line Ysteam,sat as a function of temperature (for known mass fractions of gas mixture components far from the droplet),(27)Ysteam,sat⁡≈Ysteam,sat⁡(Tsat⁡,{Ym,m=1,6¯})=Ysteam,sat⁡(Tdrop,{Ym,m=1,6¯}),
the above equation for ωH2O will be written as (28)ωH2O=2nπddropρDsteam×ln⁡[1-Ysteam,∞1-Ysteam,sat⁡(Tdrop,{Ym,m=1,6¯})].
In order to estimate Ysteam,sat for the complex (Tsat,{Ym,m=1,6¯}), one can employ the following algorithm: based on the known value of  Tsat, we find pressure and density of steam on the saturation line; mass fractions of the other components are taken in mutual proportions corresponding to the gas mixture far from the droplet; densities of the other mixture components are chosen in such a way that the pressure of the target mixture should be equal to the gas mixture pressure far from the droplet. The resulting densities of the mixture components are used to calculate their mass fractions.Note that the process of droplet falling is observed to include two modes. The first mode of droplet motion begins, when droplets issue from the nozzles of the spray system. It is characterized by initial bulk condensation of steam on the droplets with increase in their diameter and temperature. The second mode of motion begins, when the droplet temperature reaches the saturation point for the current steam pressure. This mode is characterized by droplet evaporation and decrease in droplet diameter.During the second mode, steam mass fraction is observed to grow, while the gas mixture temperature in the region, where evaporating droplets move, decreases. Such a growth of steam content in some cases may lead to secondary spontaneous bulk condensation. This process, however, will be ignored for the reasons described below.Bulk condensation requires, in particular, that the gas mixture temperature should decrease to the saturation point (and below). If the spray system is out of operation, the temperature decreases rather slowly. For this reason, even if secondary bulk condensation occurs, this process will be slow (at least, compared to the primary condensation on the falling cold droplets of water). In addition, the group of droplets considered here will be followed by another group of droplets, and so forth. The subsequent groups of droplets will first pass through the already cooled layers of the gas mixture almost without heating. Consequently, in the mixture layers, where the first group evaporates, primary condensation of steam can occur on the second group. The process of droplet flight from the time it issues from the nozzle to its complete evaporation is rather fast (at least, compared to the variation of mixture temperature, when the spray system is out of operation). For this reason, secondary bulk condensation is ignored here.Now, let us develop the model of droplet flight in a flow of gas mixture. We introduce the parameterV⃗drop, which characterizes the average velocity of droplets in a small volume of mixture. The momentum equation can be written in the following form (see also (8)):(29)∂∂t∭VρV⃗dV=∬V(ρg⃗-∇⃗p+∇⃗τ)dV,
where p is static pressure of the gas mixture. This equation is integrated over the droplet volume. We assume here that liquid density is constant, and pressure variation in the droplet volume can be ignored: (30)d(mdropV⃗drop)dt=mdropg⃗+∬Sdropτ⃗ndS.Droplet drag is estimated here using a function that approximates experimental data [14]:(31)∬Sdropτ⃗ndS=-0.5CxρSMSV⃗drop|V⃗drop|,
where Cx is the drag coefficient of the droplet; ρ is density of the gas mixture around the droplet; SMS=πddrop2/4 is midsection area. Here, just as in (30), τ⃗n is the projection of τ to the outward normal to the droplet surface. Tabulated reference parameters of the function Cx=Cx(u) (where u is the velocity of droplet motion relative to the gas environment) have been approximated by Kobyakov in [14]: (32)Cx(u)=-3.301⋅10-7[uddropν]+0.466+45.822[uddropν]-1-28.262[uddropν]-2+11.995[uddropν]-3,
where ν is the kinematic viscosity of the gas mixture. In the first approximation, we assume that u=|V⃗drop|.Thus, the equation of droplet motion can be written in the following way:(33)∂(mdropV⃗drop)∂t=mdropg⃗-πCxρddrop2V⃗drop|V⃗drop|8.
According to Newton’s third law, as the droplet moves across the gas mixture, this mixture undergoes forces opposite to the drag forces acting on the droplet. Let us write a corresponding summand for the momentum equation of the gas mixture.We consider a small volume of the gas mixture,V0. It contains (nV0) droplets. Forces acting from them are given by the following summand:(34)(-∬Sdropτ⃗ndS)(nV0)=nV0πCxρddrop2V⃗drop|V⃗drop|8.
Dividing the summand by V0 and letting V0 goes to zero yield a final expression for the summand in the momentum equation of the gas mixture describing the influence of droplets on this gas mixture:nπCxρddrop2V⃗drop|V⃗drop|/8. To establish the final form of the momentum equation for the gas mixture, one needs to add the resulting summand to the right-hand side of (8*):(8)∂(ρV⃗)∂t+∇⃗⋅(ρV⃗V⃗)=(ρ-ρatm)g⃗-∇⃗P+∇⃗⋅(μ+ρνTμτ)-23∇⃗(ρK)+nπCxρddrop2V⃗drop|V⃗drop|8.
Let us now proceed to transformations of energy equation (9*). As before, let us consider a small volume of the gas mixture V0. The power of forces experienced by the gas mixture as a result of droplet motion equals(35)∬Sdrop(τ⋅n⃗mix)⋅V⃗dropdS=-∬Sdrop(τ⋅n⃗drop)⋅V⃗dropdS,
where n⃗mix is unit vector normal to the droplet surface external to the gas mixture; n⃗drop is unit vector normal to the droplet surface external to the droplet. Let τ⃗ndrop=τ·n⃗drop. We introduce the notion of droplet surface average viscous stress: τ⃗n,averdrop=Sdrop-1∬Sdropτ⃗ndropdS. Then,(36)∬Sdrop(τ⋅n⃗mix)⋅V⃗dropdS=-∬Sdropτ⃗ndrop⋅V⃗dropdS≈-τ⃗n,averdrop⋅V⃗dropSdrop=-V⃗drop⋅∬Sdropτ⃗ndropdS=πCxρddrop2|V⃗drop|38.
The small volume V0 contains (nV0) droplets; consequently, the summand characterizing the forces exerted by the droplets in the energy equation for the volume V0 will be given by the expression: nV0πCxρddrop2|V⃗drop|3/8. We divide the summand by V0 and let it go to zero: nπCxρddrop2|V⃗drop|3/8. To write the final form of the energy equation, one needs to add the resulting summand to the right-hand side of (9*): (9)∂(ρH)∂t+∇⃗⋅(ρHV⃗)=∂P∂t+∂Q∂t+0.5Qphase_transition[ωH2O-|ωH2O|]+(Cp)H2OTsat⁡ωH2O+ρg⃗⋅V⃗+∇⃗⋅[(λ+λT)∇⃗T+μ+ρνTμτ⋅V⃗-23ρKV⃗]+∑m=1N∇⃗⋅[(Cp)mTμ+ρνTScm∇⃗Ym]+nπCxρddrop2|V⃗drop|38.The velocity of droplet motion isV⃗drop. Consequently, it can be assumed that the variable n (specific (per unit volume of gas mixture) number of droplets) also moves at a velocity of V⃗drop. Based on these considerations, one can write an equation for n: (37)∂n∂t+V⃗drop⋅∇⃗n=0.Thus, the mathematical model of spray system operation will combine (1)–(15), (20), (21), (24), (28), (32), (33), (37)).


## 4. Accounting for Film Condensation of Steam on Containment Walls

Film condensation of steam on containment walls will be taken into account here by defining corresponding boundary conditions. Note that a simplified variant of such an approach has been presented in [18].As we know, quite an acceptable estimate of the rate of heat and mass exchange for film condensation of single-component resting steam on a vertical flat wall is provided by the composite formula of Nusselt [16, 19–22]:
(38)jb=qwallQphase_transition,qwall(x)=α(x)[Tsat⁡-Twall],q̅wall(H̃)=α̅(H̃)[Tsat⁡-Twall],α(x)=[λliquid3ρliquid2gQphase_transition4μliquid(Tsat⁡-Twall)x]0.25,α̅(H̃)=43[λliquid3ρliquid2gQphase_transition4μliquid(Tsat⁡-Twall)H̃]0.25,δ(x)=[4λliquidμliquid(Tsat⁡-Twall)xQphase_transitionρliquid2g]0.25,w̅(x)=ρliquidg3μliquidδ2(x),
where jb is the mass flux of condensing steam; qwall is the heat flux during steam condensation; x is the distance from an arbitrary point on the wall to the upper film surface (vertically); H̃ is the film height; qwall(x) and α(x) are local (for the coordinate x) values of the heat flux and heat transfer coefficient; q̅wall(H̃) and α̅(H̃) are film average (for the film height H̃) values of the heat flux density and heat transfer coefficient; λliquid, ρliquid, and μliquid are the thermal conductivity, density, and dynamic viscosity of condensate (as a rule, these are considered for the temperature 0.5(Tsat+Twall)); Twall is the wall temperature; δ(x) is the local thickness of the condensate film; w̅(x) is the velocity of condensate averaged through the film thickness. The wall temperature, saturation point, and liquid properties are assumed constant. Laminar flow of the film is observed subject to the condition  Re<Recrit=400 [19], where Re is the (film) Reynolds numbers: Re=ρliquidw̅δ/μliquid.Effects of Steam Flow Velocity
Composite formula (38) has limited application, because it has been developed for fixed dry saturated steam. However, in accordance with [16], this formula remains valid for small velocities of steam flow, |V⃗|≤10m/s. The average velocity of the gas-steam mixture in the case of interest is 1 m/s. Therefore, formula (38) can be considered applicable to our model, if steam is dry and saturated [16].Effects of Steam Superheating
In case of condensation of superheated steam, one should account for the superheating heatqsteam=Csteam(Tsteam-Tsat) and replace evaporation heat Qphase_transition in the formula  by heat of superheating Qphase_transition′=Qphase_transition+qsteam, where Tsteam is the temperature of superheated steam, Csteam is heat capacity of superheated steam [16]. Temperature difference is again assumed to be equal to ΔT=Tsat-Twall.
If the wall temperature is below the saturation point, condensation of superheated steam proceeds in a way similar to the case of saturated steam. Of course, this does not mean that superheated steam becomes immediately saturated throughout the volume. Steam becomes saturated only at the wall as it cools down, while far from the wall it can remain (and remains) superheated. SinceQphase_transition′>Qphase_transition, it follows from the above formulas that heat transfer during condensation of superheated steam is somewhat higher than in case of saturated steam. During condensation of saturated dry steam, the heat flux and steam mass flux can be found using formulas (see (38)) qwall=α(Tsat-Twall), jb=qwall/Qphase_transition.
If steam is superheated, additional heat flux,qsteam_add=qsteamjb=qsteamqwall/Qphase_transition, will be withdrawn from the gas mixture. Thus, the total heat flux toward the wall will be
(39)qwall,Σ=qwall+qsteam_add=qwall(1+qsteamQphase_transition)=qwall[1+Csteam(T-Tsat⁡)Qphase_transition]
or
(40)qwall,Σ=α(Tsat⁡-Twall)[1+Csteam(T-Tsat⁡)Qphase_transition]=αΣ(T-Twall),
where αΣ is the wall-to-superheated-steam heat transfer coefficient; note that here Tsteam=T. Hence,
(41)αΣ=α[1+Csteam(T-Tsat⁡)Qphase_transition]Tsat⁡-TwallT-Twall.Effects of Noncondensing Gas Content in Steam
The presence of air and other non-condensing gases in steam significantly decreases heat transfer during condensation [16]. The reason is that only steam condenses on the cold wall, while air does not. If there is no convection, air will build up near the wall with time and significantly hamper the flow of steam to the wall.
In fact, according to Dalton’s law, the total pressurep0 of the steam-air mixture is the sum of partial pressures psteam and pair of steam and air, respectively, that is, p0=psteam+pair. As a result of steam condensation, the pressure psteam near the wall is lower than in the rest of the volume. Therefore, the pressure psteam decreases continuously as steam approaches the wall, and the closer it comes to the wall, the faster it decreases, whereas the pressure pair increases [16]. This produces high air mass fraction near the wall, which is identified as a layer that can be penetrated by steam molecules only by way of diffusion.
As a result of these processes, the temperature differenceΔT=Tsat-Twall also decreases, because the temperature of the mixture is always equal to the temperature of steam saturation at partial pressure psteam. Since psteam<p0, the temperature Tsat is always lower than the saturation point at pressure p0.
The experimental curve of the relative heat transfer coefficient as a function of relative air content in steam is provided in [16] (Figure 1). The value of γair/γsteam is plotted here in percent along the x-axis, and the ratio αair/α is plotted along the y-axis, where γair is specific weight of air, γsteam is specific weight of steam, αair is the heat transfer coefficient in the presence of air in steam, and α is the heat transfer coefficient for condensation of pure steam. As it shown in Figure 1, the heat transfer coefficient of steam containing only 1 percent of air decreases by as mush as 60 percent.
In our case of hydrogen-steam-air flow inside the containment during severe accidents, air, which consists primarily of oxygen and nitrogen, and hydrogen can be treated as a non-condensing gases. The relative content of hydrogen (even in case of its accidental release) is rather small. For this reason, the plot (see Figure1) can be considered suitable for modeling the severe accident of interest. As seen from Figure 1, the heat transfer coefficient first decreases rapidly with air content in steam, and then, starting from around γair/γsteam≈5%, this curve reaches the asymptote αair/α≈0.16. In the case of interest, the quantity γair/γsteam is much greater than 5%. As a consequence of this, we suppose that the heat transfer coefficient αnot_condensed can be estimated in the first approximation using Nusselt’s formula (38) multiplied by 0.16, that is, αnot_condensed=0.16α.
LetrH2O be the volume fraction of steam in the mixture. Then, the partial pressure of steam pH2O can be estimated as pH2O=rH2Op. The saturation point of steam can be calculated in this case as Tsat=Tsat(pH2O)=Tsat(rH2Op), where Tsat=Tsat(p) is the temperature as a function of pressure on the saturation line of water.
For numerical simulations of steam condensation on the inner walls of a reactor building, we introduce a parameterα̅hot to characterize the average heat transfer coefficient of single-component superheated resting steam on a vertical flat wall (see (38)):
(42)α̅hot(H̃,p,T)={43λliquid3ρliquid2gQphase_transition4μliquid[Tsat⁡(p)-Twall]H̃4×(1+Csteam[T-Tsat⁡(p)]Qphase_transition)Tsat⁡(p)-TwallT-Twall,ifTsat⁡(p)>Twall;0,ifTsat⁡(p)≤Twall.
The condition in (42) is introduced to constrain the range of application for the steam condensation formula. In particular, for Tsat<Twall, condensation is replaced by evaporation. In the case of interest, the process of water film evaporation from the inner walls of the reactor building is ignored. The heat transfer coefficient during steam condensation in our case can be calculated using the formula
(43)α̃={0.16α̅hot(H̃,rH2Op,T),if0.16α̅hot(H̃,rH2Op,T)>Chot;Chot,if0.16α̅hot(H̃,rH2Op,T)≤Chot,
where Chot≈5.9W/(m2K). The value of the heat transfer coefficient due to free convection (without condensation processes) equals α=Chot [16]. The condition in (43) serves to ensure that numerical heat exchange between the gas mixture and the wall does not fall below the level of heat exchange without steam condensation.The average heat flux to the wall will be equal to q̃wall=α̃(T-Twall). As part of this heat flux is spent on steam condensation (the remaining part of the heat flux is necessary for steam cooling from T to Tsat), the mass flux of condensation can be calculated using the formula
(44)j̃b=q̃wallQphase_transition[1+Csteam(T-Tsat⁡)Qphase_transition]-1.
Simulation of steam condensation on containment walls involves solving conjugate problems of gas mixture distribution in the containment volume (1)–(15), (20), (21), (24), (28), (32), (33), (37)] and heat conduction in the wall:
(45)ρwallCwall∂Twall∂t-∇⃗⋅(λwall∇⃗Twall)=0,
where ρwall, Cwall, Twall, and λwall are density, specific (per unit mass) heat capacity, temperature, and thermal conductivity of the wall. For the wall/condensate interface one should assign the Neumann boundary conditions:
(46)λwall∂T∂n=q̃wall.
Here and below, n⃗ is unit outward normal to the surface of the region of interest. On the gas mixture side of the same interface one should use the Neumann boundary conditions for steam “withdrawal”:
(47)ρDsteam∂(Yatm_H2O+YH2O+YPAR_prod)∂n=-j̃b.
By analogy with (47), the boundary condition for the energy equation will be
(48)λ∂T∂n=-q̃wall,
to characterize the heat flux generated by gas-to-liquid “transition” of water (for the mass flux j̃b) at Tsat.
The saturation point of steam can be calculated using the equation similar to the Antoine equation [23]:
(49)p(Tsat⁡)=133.32exp⁡[18.76-4081.18Tsat⁡-36.38]
or
(50)Tsat⁡(p)=36.38-4081.18[ln⁡(p133.32)-18.76]-1,
where dimensionality in the formula has the form of [p]=[Pa], [Tsat]=[K].Figure 1
Relative variation of the heat transfer coefficient during condensation as a function of air content in steam (subscript “b” stands for “air”; subscript “Π” stands for “steam”) [16].


## 5. Accounting for the Effects of the Containment’s Steel Liner

In new NPP designs, the inside surface of double containment inner concrete wall has a liner up to 10 mm thick carbon steel plates. The thickness of inner wall can vary from one meter to a meter and a half. In order to determine the appropriate mesh density for the numerical analysis of the set of (1)–(15), (20), (21), (24), (28), (32), (33), (37), one should first estimate conventional contributions of each of stated structural parts  to the process of heat withdrawal. Let us present one possible way for estimating such conventional contributions.The first step in this method is to perform a preliminary numerical analysis of model (1)–(15), (20), (21), (24), (28), (32), (33), (37) with adiabatic boundary conditions on containment walls. It should be noted here that this set of equations completed by relevant boundary conditions can be analyzed successfully using the CFD code ANSYS/CFX [24]. In this case, the set of (1)–(15), (20), (21), (24), (28), (32), (33), (37) will need to be slightly modified without loss of originally achieved completeness and adequacy of fluid dynamics modeling of processes inside the containment during severe accidents. Such modifications are technical, and their description is therefore beyond the scope of this paper.The mentioned code has been chosen because the CFX solver [24], which is part of ANSYS/CFX, has demonstrated satisfactory results during its verification by test simulations of the standard international problem ISP-47 and its analogs under an IAEA task order  [2, 13, 25].The second step is solving a model heat transfer problem, the sketch of which is shown in Figure2.Figure 2
Sketch of a model heat transfer problem to estimate the parameters of heat exchange through the inner containment wall.This problem models a transient heat flow through some region of a nonuniform wall composed of steel and concrete layers of corresponding thicknesses. Considering the size of the reactor building structure, the curvature of the containment’s cylindrical and spherical parts can be ignored when making the needed estimates, and the problem can be solved in the plane statement. The heat transfer problem can be solved by the Finite Element Method (FEM) in the program ANSYS environment [26].In the numerical thermal analysis of the model problem the Newton boundary conditions were assigned on the outside (concrete) and inside (steel) surfaces of the wall. Let us give an example of values of the physical parameters we used for the numerical analysis of an severe accident at the Baltic NPP. The following conditions were defined on the outside wall surface: ambient temperatureTam=293K; heat transfer coefficient for free convection αout=Chot. On the wall inside surface the gas mixture temperature Tmix(t) varied by the law, depicted in Figure 3. For the heat transfer coefficient was taken equality αin=αout. Such a value of αin was chosen based on the results of preliminary CFD analysis of steam-hydrogen mixture flow in reactor building rooms (Figure 4).Figure 3
Volume-averaged temperature of the steam-hydrogen-air mixture for adiabatic boundary conditions.Figure 4
Example of a snapshot of the field of velocities (m/s) in containment rooms.Physical properties of the materials were specified in accordance with [27]: steel and concrete densities ρsteel=7850kg/m3 and ρconcrete=2400kg/m3specific heat capacities of steel and concrete Csteel=460J/(kg·K) and Cconcrete=780J/(kg·K); thermal conductivities of steel and concrete λsteel=45W/(m·K) and λconcrete=1.23W/(m·K). The difference in physical properties of reinforcement, embedded parts, and other “nonconcrete” structural members of the containment concrete wall was ignored, because the volume of these elements is small compared to the volume of concrete.Cross-sections of the computational model (see Figure2) were adiabatic to represent symmetry about boundaries of an arbitrarily chosen part of the extended wall. The model shown in Figure 2 can be of any height, because comparison is performed for specific thermal characteristics  of the containment’s structural parts.These numerical simulations produced time distributions of temperature through the thickness of the containment wall. Figure5 shows the through-the-thickness temperature distribution in the inner containment wall at the time 10 000 s after the beginning of an accidental steam-hydrogen mixture release. As we see in Figure 5, only a small part of the wall structure can heat up over the time interval of interest. The depth of heat penetration into concrete can be quantified using Figure 6, which shows through-the-thickness temperature field profiles for several time points.Figure 5
Temperature pattern (K) through the thickness of the inner containment wall at the time 10,000 s after the beginning of accidental steam-hydrogen mixture release.Figure 6
Temporal variations of the through-the-thickness temperature profile of the inner containment wall during accidental steam-hydrogen mixture release.It follows from Figure6 that the maximum depth of heat penetration in the wall during the time interval of 10,000 s from the beginning of the severe accident does not exceed 250÷300mm. In the wall’s concrete layers situated far from the heated surface, temperature variation during the accident stays within ΔT≪1K. Therefore, for the numerical analysis of severe accidents in the containment wall/environment heat exchange it would suffice to supplement the computational model with solid domain consists of steel liner and 250 mm concrete layer. On the outside surface of the concrete layer one should assign the Dirichlet boundary condition: T=Tam.The simulation results shown in Figure5 are evidence of that the maximum temperature arises in the thin layer of the wall’s steel liner. At the same time, although the temperatures in concrete are much lower, the depth of heat penetration in it is much greater than the thickness of steel liner. In order to evaluate the possibility of ignoring the contribution of the steel liner to heat withdrawal (to optimize the model size), one should compare these heat quantities accumulated with time in the concrete layer and steel liner. Heat accumulated on unit area of the side surface of the containment wall by the time t equals(51)Q(t)=∫x1x2ρC[T(x,t)-T(x,t0)]dx,
where T(x,t) is the function of through-the-thickness temperature distribution in the region [x1;x2] with uniform physical properties ρ and C; t0 is initial time.Numerical integrating of the results of solving the model problem for each of the containment wall materials in accordance with the above formula for a number of timest∈(0;10000s] for the problem-specific fixed values of parameters (t0=0; x1=0, x2=0.006m for steel liner; x1=0.006m, x2=0.3m for concrete wall; ρ and C corresponding to the above values for steel and concrete) gives time profiles of heat accumulated in the wall’s structural parts. These profiles are shown in Figure 7.Figure 7
Heat accumulated in the wall’s structural parts during the severe accident.As evidenced by the numerical simulation results (see Figure7), the basic contribution to the process of heat withdrawal from the gas mixture inside the containment is made by the containment’s concrete wall. However, the fraction of heat accumulated in steel liner is also considerable. The minimum value of this fraction is observed at the end of the accident time interval (t=10000s) and equals 17 percent of concrete heat. Thus, it is inadmissible to ignore the contribution of the containment wall’s steel liner to heat exchange with the environment in high accurate simulations.On the other hand, for direct numerical simulations of steel liner, the model should include a fine mesh with a characteristic linear dimension ofle≤0.006m extending over the entire surface area of the containment. As a result, the mesh density for the problem of interest will increase unreasonably. In addition, in accordance with the known Courant criterion, the maximum time step should also be reduced substantially in order to provide convergence of the numerical solution of the transient problem. Thus, it does not seem possible to use direct numerical simulations of the wall’s steel liner for the numerical analysis of the severe accident of interest under the current limitations on computer and time resources available in design organizations.At the same time, because of steel’s high thermal conductivity, the temperature difference between the outer (containment facing) and inner (concrete facing) surfaces of the liner structure is extremely small (see Figures5 and 6). This difference can therefore be ignored, and steel liner can be assumed to have uniform through-the-thickness temperature at each point of time (model of thermal conductive thin shell). This assumption alone makes it possible to simulate the two-layer containment wall (see Figure 2) without remeshing the model and reducing the time step in most of commercial CFD codes.Examples of CFD analysis of steam-hydrogen-air cloud evolution inside the Baltic NPP containment with a functioning spray system, PARs, and condensation heat exchangers of passive heat removal from the containment rooms are illustrated by some results shown in Figure8.Calculated results of the temperature field (K) (a) and isosurface of 40-percent molar fraction of steam (b) in containment rooms.
(a)(b)The modeling area can be described by the following principal dimension: the height is about 70 m; the radius of a spherical segment is 22 m; the height of a parallel portion is about 18 m; the radius of a parallel portion is 22 m. Here were analyzed the consequences from hypothetical accident due to small leak of the equivalent diameter 100 mm from the cold line of the coolant pipe at 6.9 m distance from the reactor inlet and with the superposition of active part fault of the area emergency cooling system. The inflow functions of hydrogen and steam to the modeling area are presented in Figure9.Integral operation (time) of  hydrogen (a) and steam (b).
(a)(b)


## 6. About Models Validation

The examples of 3D CFD analysis results are presented in Figures4 and 8. The authors donot dispose of the information about the implementation of full-scale experiments in that problem formulation somewhere in the world. Taking into account a big volume of the containment (height 70 m, base diameter 44 m) and also the mass of realised steam-hydrogen mixture (up to 120 tons) the implementation of full-scale experiments, such as severe accident inside the reactor building of NPP, seems to be very difficult. That is why mathematical modeling is widespread for the analysis of severe accidents at NNP nowadays. The adequacy of mathematical models and methods for numerical simulations is proved in theory and by the results of specially formulated laboratory experiments. Such approach is particularly described in monograph [28].It should be noted that the relevance of(k-ω) turbulence model in (10, 11) with regard to steam flow was done in comparison of the results of numerical analysis with data from [29–32]. The comparison of computational and experimental results was done for such parameters, as components of flow velocity, concentrations of mixture species, density, and temperature. Validation and calibration of PAR computational models were carried out on the basis of published data, obtained by empirical correlations (hydrogen consumption rate), and from the laboratory experiment (temperature, composition of gas mixtures) and CFD-analysis [3–5, 33]. The verification of the above stated spray system models was done in the same way. For the comprehensive verification of proposed models of spray systems were used a wide range of theoretical, experimental and calculation data, which were presented in [1–3, 6, 12, 13, 29–32].The description of numerical experiments and comparison procedure is beyond the scope of this paper. However, it should be noted that the verification outcome has shown the reasonability of proposed models application for CFD-analysis of the emergency emissions of steam-hydrogen-air mixture inside containment.


## 7. Conclusions

The presented in the paper model of steam-hydrogen-air  mixture evolution in containment rooms accounting for the operation of the spray system enables high accurate numerical analysis of severe accidents without high-performance computers. With properly specified boundary conditions, this model makes it possible to perform numerical analysis accounting for the operation of PARs and condensation heat exchangers of the passive heat removal system in the reactor building. It is reasonable to use this model as a substitute for simplified models in alike simulations described in [1, 18].



---

*Source: 284274-2011-11-14.xml*