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.

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, (𝑚=1); air as a mixture of its three basic components—oxygen, (𝑚=2), nitrogen, (𝑚=3), and atmospheric water vapor, (𝑚=4); steam injected into the containment during severe accidents, (𝑚=5); additional steam produced by hydrogen and oxygen absorption by PAR [4–6] (referred to as absorption product) (𝑚=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 (𝑘−𝜔)-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 (𝑘−𝜀)-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]):𝜕𝜌+⃗𝜌⃗𝐕𝜕𝑡∇⋅=0,(1∗)

𝜕(𝜌𝐾)+⃗⃗𝐕=⃗𝜕𝑡∇⋅𝜌𝐾∇⋅𝜇+𝜌𝜈𝑇Pr𝐾⃗∇𝐾+𝐺−𝛽∗𝜈𝜌𝐾𝜔−𝑇Pr𝑇⃗⃗𝐠⋅∇𝜌,(10)𝜕(𝜌𝜔)+⃗⃗𝐕=⃗𝜕𝑡∇⋅𝜌𝜔∇⋅𝜇+𝜌𝜈𝑇Pr𝜔⃗+𝜔∇𝜔𝐾𝛼𝐺−𝛽∗∗𝜈𝜌𝐾𝜔+𝑇Pr𝑇×⃗⃗−⃗⃗,(𝛼+1)max𝐠⋅∇𝜌;0𝐠⋅∇𝜌(11)𝑃=𝜌𝑅0𝑇𝑁𝑚=1𝑌𝑚𝑀𝑚+𝜌atm𝑔𝑥3−𝑥3,0,𝐻=𝐶𝑝⃗⃗𝑇+0.5𝐕⋅𝐕;𝛼𝑆=𝑁𝑚=1𝑓𝛼𝑌𝑚,𝛼𝑚,𝛼≡𝜇,𝜆,𝐶𝑝,𝑇(12)𝜇=273.153/2273.15+𝐶𝑆𝑇+𝐶𝑆𝜇0,𝜆𝑇=𝐶𝑝𝜌𝜈𝑇Pr𝑇;Sc𝑚=𝜇𝜌𝐷𝑚,𝜏(13)𝑖𝑗=𝜇𝜕𝑢𝑖𝜕𝑥𝑗+𝜕𝑢𝑗𝜕𝑥𝑖−23𝛿𝑖𝑗𝜕𝑢𝑘𝜕𝑥𝑘,(14)𝐺=𝜌𝜈𝑇12𝜕𝑢𝑖𝜕𝑥𝑗+𝜕𝑢𝑗𝜕𝑥𝑖2−23𝜕𝑢𝑘𝜕𝑥𝑘2−23𝜌𝐾𝜕𝑢𝑘𝜕𝑥𝑘,𝜈𝑇=𝛼1𝐾𝛼max1𝜔;𝑆𝐹2,𝑆=2𝑆𝑖𝑗𝑆𝑖𝑗,𝑆𝑖𝑗=12𝜕𝑢𝑖𝜕𝑥𝑗+𝜕𝑢𝑗𝜕𝑥𝑖,𝐹2=tanhmax22√𝐾𝛽∗;𝜔𝑦500𝑣𝑦2𝜔,(15)
where ⃗⃗𝐷(⋯)/𝐷𝑡≡(𝜕(⋯)/𝜕𝑡)+𝐕⋅∇(⋯) is the substantial derivative of the scalar function (writing a substantial derivative of a vector function means substantial differentiation of vector function components); 𝑡 is time; ⃗𝐕 is velocity with components 𝑢1, 𝑢2, 𝑢3; ⃗∇ is nabla; 𝜌, 𝜌𝑚, 𝜌atm are density of mixture, 𝑚th mixture component, and atmosphere at some point in space of containment rooms, respectively; 𝑌H2=𝑌1, 𝑌O2=𝑌2, 𝑌N2=𝑌3, 𝑌atm_H2O=𝑌4, 𝑌H2O=𝑌5, 𝑌PAR_prod=𝑌6, 𝑌𝑚=𝜌𝑚/𝜌 is the relative mass fraction of the 𝑚th component; 𝜇 is the dynamic viscosity; 𝜇𝑇=𝜌𝜈𝑇 is the turbulent viscosity; 𝜈 is the kinematic viscosity; Sc𝑚 is the Schmidt number corresponding to the 𝑚th component; 𝜔H2 is the mass rate of the generalized irreversible single-stage reaction of hydrogen absorption in PAR [4, 5]; 𝜉Omass2 is the stoichiometric mass coefficient of oxygen; 𝑁 is the number of components in the gas mixture (in our case, 𝑁=6 (see above)); the subscript 𝑇 stands for “turbulent,” and the subscript atm, for the “atmosphere in containment rooms”; ⃗𝐠 is gravity acceleration (𝑔 is the gravity acceleration modulus); 𝑃 is piezometric pressure of mixture; 𝝉 is tensor of viscous stresses with components 𝜏𝑖𝑗; 𝐾 is turbulent kinetic energy; ⃗⃗𝐕𝐻=ℎ+0.5𝐕⋅ is total specific (per unit mass) enthalpy (ℎ=𝐶𝑝𝑇 is enthalpy for perfect gas, 𝐶𝑝 is specific (per unit mass) heat capacity at constant pressure; 𝑇 is mixture temperature); 𝜕𝑄/𝜕𝑡 is the specified rate of heat release from outer sources per unit volume (heat generation rate per unit volume); 𝜆 is the thermal conductivity; 𝜆𝑇 is the turbulent thermal conductivity; Pr is the Prandtl number (the subscript 𝐾 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)); Pr𝐾=2.0 [7–11]; 𝛽∗=0.09 [7–9]; 𝜔 are turbulence frequencies; Pr𝜔=2 [7–9]; 𝛼=5/9 [7–9]; 𝛽∗∗=0.075 [7–9]; 𝑅0 is the universal gas constant; 𝑀𝑚 is molar mass of the 𝑚th component; 𝑥3 is the vertical coordinate in the Cartesian frame of reference directed away from the earth center (here, 𝑥1, 𝑥2, 𝑥3 are radius vector coordinates of the point in the Cartesian frame of reference); 𝑥3,0 is the coordinate of the point 𝑥3 with a defined value of 𝜌atm; 𝑓𝛼(⋯) are known semiempirical functions; 𝐶𝑆 is the Sutherland constant; 𝜇0 is the dynamic viscosity under normal conditions; 𝐷𝑚 is the coefficient of binary diffusion of the mth component in the rest of the mixture; 𝛿𝑖𝑗 is the Kronecker delta; 𝛼1=5/9 [7–9]; 𝑦 is the distance to the nearest wall. The functions 𝐺, 𝑆, 𝐹2 in (15) are auxiliary. Note that the model uses Favre-averaged velocity components and thermal variables (𝐻, ℎ, 𝑇) [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: 𝜕𝜌+⃗𝜌⃗𝐕𝜕𝑡∇⋅=𝜔H2O;(1)

𝜕𝜌𝑌H2O+⃗𝜕𝑡∇⋅𝜌𝑌H2O⃗𝐕=⃗∇⋅𝜇+𝜌𝜈𝑇ScH2O⃗∇𝑌H2O+𝜔H2O;(6)

𝜕(𝜌𝐻)+⃗⃗𝐕=𝜕𝑡∇⋅𝜌𝐻𝜕𝑃+𝜕𝑡𝜕𝑄𝜕𝑡+0.5𝑄phase_transition×𝜔H2O−||𝜔H2O||+𝐶𝑝H2O𝑇sat𝜔H2O⃗⃗𝐕+⃗+𝜌𝐠⋅∇⋅𝜆+𝜆𝑇⃗∇𝑇+𝜇+𝜌𝜈𝑇𝜇⃗2𝝉⋅𝐕−3⃗𝐕+𝜌𝐾𝑁∑𝑚=1⃗𝐶∇⋅𝑝𝑚𝑇𝜇+𝜌𝜈𝑇Sc𝑚⃗∇𝑌𝑚,(9∗∗)
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; 𝑄phase_transition is specified latent mass heat of steam condensation into water droplets at 𝜔H2O<0; 𝑇sat 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 evaporation 𝑄phase_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.5𝑄phase_transition[𝜔H2O−|𝜔H2O|]. The summand in (9**) given by [(𝐶𝑝)H2O𝑇sat𝜔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 𝐽 from a spherical droplet: 𝐽=𝐽1=𝜌𝐷steam0.5𝑑dropln1−𝑌steam,∞1−𝑌steam,sat,(16)𝐽=𝐽2=𝜆steam0.5𝑑drop𝐶𝑝H2O𝐶ln1+𝑝H2O𝑇∞−𝑇sat𝑄phase_transition,(17)
where 𝜌 is gas mixture density; 𝐷steam is the coefficient of steam diffusion in the mixture; 𝑑drop is droplet diameter; 𝑌steam,∞ and 𝑇∞ are the mass fraction and temperature of steam far from the surface, respectively (where effects produced by the evaporating steam can be ignored); 𝑌steam,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 𝑆drop=𝜋𝑑2drop, 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:𝑑𝑚steam_drop𝑑𝑡=−𝑑𝑚drop𝑑𝑡=𝜋𝑑2drop𝐽,(18)
where 𝑑𝑚steam_drop/𝑑𝑡 is the rate of steam production from one droplet; 𝑚drop is droplet mass; 𝑑𝑚drop/𝑑𝑡 is the rate of droplet mass variation. Let 𝑛 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)) 𝜔H2O=𝑛𝑑𝑚steam_drop𝑑𝑡=−𝑛𝑑𝑚drop𝑑𝑡=𝑛𝜋𝑑2drop𝐽=2𝑛𝜋𝑑drop𝜆steam𝐶𝑝H2O𝐶ln1+𝑝H2O𝑇−𝑇sat𝑄phase_transition.(19)
Here, as with (1)–(15), the temperature 𝑇 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 𝑚drop=𝜌liquid𝑉drop=𝜌liquid𝜋𝑑3drop6,(20)
where 𝜌liquid is density of the droplet liquid (water); 𝑉drop is droplet volume. Hence,𝑑drop=[6𝑚drop/(𝜋𝜌liquid)]1/3, and
𝑑𝑑𝑑𝑡drop=136𝜋𝜌liquid1/3𝑚−2/3drop⋅𝑑𝑚drop𝜔𝑑𝑡=−H2O63𝑛𝜋𝜌liquid1/3𝑚−2/3drop.(21)
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: 𝑇sat≈𝑇drop. 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), ℎliquid=𝜌liquid𝐶liquid𝑇liquid, where ℎliquid is specific (per unit volume) enthalpy of the liquid; 𝜌liquid=const is liquid density (the liquid is incompressible); 𝐶liquid=const is specific (per unit mass) heat capacity of the liquid; 𝑇liquid is liquid temperature. Considering this EOS, droplet enthalpy can be written as 𝐻drop=𝑚drop𝐶liquid𝑇drop=𝜌liquid𝑉drop𝐶liquid𝑇drop, where 𝐻drop and 𝑇drop are average enthalpy and temperature of the droplet, respectively. Hence, 𝑑𝐻𝑑𝑡drop=𝐶liquid𝑚drop𝑑𝑇drop𝑑𝑡+𝑇drop𝑑𝑚drop𝑑𝑑𝑡or𝐻𝑑𝑡drop=𝐶liquid𝑚drop𝑑𝑇drop𝑑𝑡−𝑇drop𝜔H2O𝑛.(22)
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 𝑇sat:𝑑𝐻𝑑𝑡drop=𝑄phase_transition+𝐶liquid𝑇sat𝑑𝑚drop𝑄𝑑𝑡=−phase_transition+𝐶liquid𝑇sat𝜔H2O𝑛.(23)
By equating the right-hand sides of the second equation in (22) and (23), we have𝑑𝑇drop=𝜔𝑑𝑡H2O𝑛𝑚drop𝑇drop−𝑇sat−𝑄phase_transition𝐶liquid𝜔=−H2O𝑄phase_transition𝑛𝑚drop𝐶liquid.(24)

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∇⋅𝜌𝑤𝑌steam−𝜌𝐷steam∇𝑌steam=0(25)
for the following boundary conditions: 𝑟=𝑟0=0.5𝑑drop∶𝑌steam=𝑌steam,sat, 𝑟=∞∶𝑌steam=𝑌steam,∞, where 𝑟 is distance from the droplet center to an arbitrary point.

During droplet evaporation, 𝑌steam,∞≤𝑌steam,sat, which produces a flow of steam from the droplet surface to the gas mixture. If the situation is opposite, 𝑌steam,∞>𝑌steam,sat (i.e., the steam mass fraction in the gas mixture is higher than on the droplet surface equal to 𝑌steam,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, 𝜔H2O=𝑛𝜋𝑑2drop𝐽=2𝑛𝜋𝑑drop𝜌𝐷steamln1−𝑌steam,∞1−𝑌steam,sat.(26)
If we first write the mass fraction of steam on the saturation line 𝑌steam,sat as a function of temperature (for known mass fractions of gas mixture components far from the droplet),𝑌steam,sat≈𝑌steam,sat𝑇sat,𝑌𝑚,𝑚=1,6=𝑌steam,sat𝑇drop,𝑌𝑚,𝑚=,1,6(27)
the above equation for 𝜔H2O will be written as 𝜔H2O=2𝑛𝜋𝑑drop𝜌𝐷steam⎡⎢⎢⎣×ln1−𝑌steam,∞1−𝑌steam,sat𝑇drop,𝑌𝑚,𝑚=⎤⎥⎥⎦.1,6(28)
In order to estimate 𝑌steam,sat for the complex (𝑇sat,{𝑌𝑚,𝑚=1,6}), one can employ the following algorithm: based on the known value of 𝑇sat, 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 parameter 𝐕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)
where 𝑝 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: 𝑑𝑚drop𝐕drop𝑑𝑡=𝑚drop⃗𝐠+𝑆drop⃗𝝉𝑛𝑑𝑆.(30)

Droplet drag is estimated here using a function that approximates experimental data [14]:𝑆drop⃗𝝉𝑛𝑑𝑆=−0.5C𝑥𝜌𝑆MS𝐕drop||𝐕drop||,(31)
where 𝐶𝑥 is the drag coefficient of the droplet; 𝜌 is density of the gas mixture around the droplet; 𝑆MS=𝜋𝑑2drop/4 is midsection area. Here, just as in (30), ⃗𝝉𝑛 is the projection of 𝝉 to the outward normal to the droplet surface. Tabulated reference parameters of the function 𝐶𝑥=𝐶𝑥(𝑢) (where 𝑢 is the velocity of droplet motion relative to the gas environment) have been approximated by Kobyakov in [14]: 𝐶𝑥(𝑢)=−3.301⋅10−7𝑢𝑑drop𝜈+0.466+45.822𝑢𝑑drop𝜈−1−28.262𝑢𝑑drop𝜈−2+11.995𝑢𝑑drop𝜈−3,(32)
where 𝜈 is the kinematic viscosity of the gas mixture. In the first approximation, we assume that 𝐕𝑢=|drop|.

Thus, the equation of droplet motion can be written in the following way:𝜕𝑚drop𝐕drop𝜕𝑡=𝑚drop⃗𝐠−𝜋𝐶𝑥𝜌𝑑2drop𝐕drop||𝐕drop||8.(33)
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, 𝑉0. It contains (𝑛𝑉0) droplets. Forces acting from them are given by the following summand:−𝑆drop⃗𝝉𝑛𝑑𝑆𝑛𝑉0=𝑛𝑉0𝜋𝐶𝑥𝜌𝑑2drop𝐕drop||𝐕drop||8.(34)
Dividing the summand by 𝑉0 and letting 𝑉0 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:𝑛𝜋𝐶𝑥𝜌𝑑2drop𝐕drop|𝐕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*):𝜕𝜌⃗𝐕+⃗𝜌⃗𝐕⃗𝐕=𝜕𝑡∇⋅𝜌−𝜌atm⃗⃗⃗𝐠−∇𝑃+∇⋅𝜇+𝜌𝜈𝑇𝜇𝝉−23⃗∇(𝜌𝐾)+𝑛𝜋𝐶𝑥𝜌𝑑2drop𝐕drop||𝐕drop||8.(8)
Let us now proceed to transformations of energy equation (9*). As before, let us consider a small volume of the gas mixture 𝑉0. The power of forces experienced by the gas mixture as a result of droplet motion equals𝑆drop𝐧𝝉⋅mix⋅𝐕drop𝑑𝑆=−𝑆drop𝐧𝝉⋅drop⋅𝐕drop𝑑𝑆,(35)
where 𝐧mix is unit vector normal to the droplet surface external to the gas mixture; 𝐧drop is unit vector normal to the droplet surface external to the droplet. Let ⃗𝝉drop𝑛𝐧=𝝉⋅drop. We introduce the notion of droplet surface average viscous stress: ⃗𝝉drop𝑛,aver=𝑆−1drop∬𝑆drop⃗𝝉drop𝑛𝑑𝑆. Then,𝑆drop𝐧𝝉⋅mix⋅𝐕drop𝑑𝑆=−𝑆drop⃗𝝉drop𝑛⋅𝐕drop𝑑𝑆≈−⃗𝝉drop𝑛,aver⋅𝐕drop𝑆drop𝐕=−drop⋅𝑆drop⃗𝝉drop𝑛=𝑑𝑆𝜋C𝑥𝜌𝑑2drop||𝐕drop||38.(36)
The small volume 𝑉0 contains (𝑛𝑉0) droplets; consequently, the summand characterizing the forces exerted by the droplets in the energy equation for the volume 𝑉0 will be given by the expression: 𝑛𝑉0𝜋𝐶𝑥𝜌𝑑2drop|𝐕drop|3/8. We divide the summand by 𝑉0 and let it go to zero: 𝑛𝜋𝐶𝑥𝜌𝑑2drop|𝐕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*): 𝜕(𝜌𝐻)+⃗⃗𝐕=𝜕𝑡∇⋅𝜌𝐻𝜕𝑃+𝜕𝑡𝜕𝑄𝜕𝑡+0.5𝑄phase_transition𝜔𝐻2𝑂−||𝜔𝐻2𝑂||+𝐶𝑝H2O𝑇sat𝜔H2O⃗⃗𝐕+⃗+𝜌𝐠⋅∇⋅𝜆+𝜆𝑇⃗∇𝑇+𝜇+𝜌𝜈𝑇𝜇⃗2𝝉⋅𝐕−3⃗𝐕+𝜌𝐾𝑁𝑚=1⃗𝐶∇⋅𝑝𝑚𝑇𝜇+𝜌𝜈𝑇Sc𝑚⃗∇𝑌𝑚+𝑛𝜋C𝑥𝜌𝑑2drop||𝐕drop||38.(9)

The velocity of droplet motion is 𝐕drop. Consequently, it can be assumed that the variable 𝑛 (specific (per unit volume of gas mixture) number of droplets) also moves at a velocity of 𝐕drop. Based on these considerations, one can write an equation for 𝑛: 𝜕𝑛+𝐕𝜕𝑡drop⋅⃗∇𝑛=0.(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]:
𝑗𝑏=𝑞wall𝑄phase_transition,𝑞wall𝑇(𝑥)=𝛼(𝑥)sat−𝑇wall,𝑞wall𝐻=𝛼𝐻𝑇sat−𝑇wall,𝜆𝛼(𝑥)=3liquid𝜌2liquid𝑔𝑄phase_transition4𝜇liquid𝑇sat−𝑇wall𝑥0.25,𝛼𝐻=43𝜆3liquid𝜌2liquid𝑔𝑄phase_transition4𝜇liquid𝑇sat−𝑇wall𝐻0.25,𝛿(𝑥)=4𝜆liquid𝜇liquid𝑇sat−𝑇wall𝑥𝑄phase_transition𝜌2liquid𝑔0.25,𝜌𝑤(𝑥)=liquid𝑔3𝜇liquid𝛿2(𝑥),(38)
where 𝑗𝑏 is the mass flux of condensing steam; 𝑞wall is the heat flux during steam condensation; 𝑥 is the distance from an arbitrary point on the wall to the upper film surface (vertically); 𝐻 is the film height; 𝑞wall(𝑥) and 𝛼(𝑥) are local (for the coordinate 𝑥) values of the heat flux and heat transfer coefficient; 𝑞wall(𝐻) and 𝛼(𝐻) are film average (for the film height 𝐻) 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(𝑇sat+𝑇wall)); 𝑇wall is the wall temperature; 𝛿(𝑥) is the local thickness of the condensate film; 𝑤(𝑥) 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=𝜌liquid𝑤𝛿/𝜇liquid.

Effects of Steam Flow VelocityComposite 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, |⃗𝐕|≤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 SuperheatingIn case of condensation of superheated steam, one should account for the superheating heat 𝑞steam=𝐶steam(𝑇steam−𝑇sat) and replace evaporation heat 𝑄phase_transition in the formula by heat of superheating 𝑄phase_transition=𝑄phase_transition+𝑞steam, where 𝑇steam is the temperature of superheated steam, 𝐶steam is heat capacity of superheated steam [16]. Temperature difference is again assumed to be equal to Δ𝑇=𝑇sat−𝑇wall. 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. Since 𝑄phase_transition>𝑄phase_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)) 𝑞wall=𝛼(𝑇sat−𝑇wall), 𝑗𝑏=𝑞wall/𝑄phase_transition.If steam is superheated, additional heat flux, 𝑞steam_add=𝑞steam𝑗𝑏=𝑞steam𝑞wall/𝑄phase_transition, will be withdrawn from the gas mixture. Thus, the total heat flux toward the wall will be
𝑞wall,Σ=𝑞wall+𝑞steam_add=𝑞wall𝑞1+steam𝑄phase_transition=𝑞wall𝐶1+steam𝑇−𝑇sat𝑄phase_transition(39)
or
𝑞wall,Σ𝑇=𝛼sat−𝑇wall𝐶1+steam𝑇−𝑇sat𝑄phase_transition=𝛼Σ𝑇−𝑇wall,(40)
where 𝛼Σ is the wall-to-superheated-steam heat transfer coefficient; note that here 𝑇steam=𝑇. Hence,
𝛼Σ𝐶=𝛼1+steam𝑇−𝑇sat𝑄phase_transition𝑇sat−𝑇wall𝑇−𝑇wall.(41)

Effects of Noncondensing Gas Content in SteamThe 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 pressure 𝑝0 of the steam-air mixture is the sum of partial pressures 𝑝steam and 𝑝air of steam and air, respectively, that is, 𝑝0=𝑝steam+𝑝air. As a result of steam condensation, the pressure 𝑝steam near the wall is lower than in the rest of the volume. Therefore, the pressure 𝑝steam decreases continuously as steam approaches the wall, and the closer it comes to the wall, the faster it decreases, whereas the pressure 𝑝air 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 Δ𝑇=𝑇sat−𝑇wall also decreases, because the temperature of the mixture is always equal to the temperature of steam saturation at partial pressure 𝑝steam. Since 𝑝steam<𝑝0, the temperature 𝑇sat is always lower than the saturation point at pressure 𝑝0.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 Figure 1) 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𝛼.Let 𝑟H2O be the volume fraction of steam in the mixture. Then, the partial pressure of steam 𝑝H2O can be estimated as 𝑝H2O=𝑟H2O𝑝. The saturation point of steam can be calculated in this case as 𝑇sat=𝑇sat(𝑝H2O)=𝑇sat(𝑟H2O𝑝), where 𝑇sat=𝑇sat(𝑝) 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)):
𝛼hot=⎧⎪⎪⎪⎨⎪⎪⎪⎩4𝐻,𝑝,𝑇34⎷𝜆3liquid𝜌2liquid𝑔𝑄phase_transition4𝜇liquid𝑇sat(𝑝)−𝑇wall𝐻×𝐶1+steam𝑇−𝑇sat(𝑝)𝑄phase_transition𝑇sat(𝑝)−𝑇wall𝑇−𝑇wall,if𝑇sat(𝑝)>𝑇wall;0,if𝑇sat(𝑝)≤𝑇wall.(42)
The condition in (42) is introduced to constrain the range of application for the steam condensation formula. In particular, for 𝑇sat<𝑇wall, 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
⎧⎪⎪⎨⎪⎪⎩𝛼=0.16𝛼hot𝐻,𝑟H2O,𝑝,𝑇if0.16𝛼hot𝐻,𝑟H2O𝑝,𝑇>𝐶hot;𝐶hot,if0.16𝛼hot𝐻,𝑟H2O𝑝,𝑇≤𝐶hot,(43)
where 𝐶hot≈5.9W/(m2K). The value of the heat transfer coefficient due to free convection (without condensation processes) equals 𝛼=𝐶hot [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 ̃𝑞wall=𝛼(𝑇−𝑇wall). As part of this heat flux is spent on steam condensation (the remaining part of the heat flux is necessary for steam cooling from 𝑇 to 𝑇sat), the mass flux of condensation can be calculated using the formula
̃𝑗𝑏=̃𝑞wall𝑄phase_transition𝐶1+steam𝑇−𝑇sat𝑄phase_transition−1.(44)
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:
𝜌wall𝐶wall𝜕𝑇wall−⃗𝜆𝜕𝑡∇⋅wall⃗∇𝑇wall=0,(45)
where 𝜌wall, 𝐶wall, 𝑇wall, 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:
𝜆wall𝜕𝑇𝜕𝑛=̃𝑞wall.(46)
Here and below, ⃗𝐧 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”:
𝜌𝐷steam𝜕𝑌atm_H2O+𝑌H2O+𝑌PAR_prod̃𝑗𝜕𝑛=−𝑏.(47)
By analogy with (47), the boundary condition for the energy equation will be
𝜆𝜕𝑇𝜕𝑛=−̃𝑞wall,(48)
to characterize the heat flux generated by gas-to-liquid “transition” of water (for the mass flux ̃𝑗𝑏) at 𝑇sat. The saturation point of steam can be calculated using the equation similar to the Antoine equation [23]:
𝑝𝑇sat=133.32exp18.76−4081.18𝑇sat−36.38(49)
or
𝑇sat𝑝(𝑝)=36.38−4081.18ln133.32−18.76−1,(50)
where dimensionality in the formula has the form of [𝑝]=[𝑃𝑎], [𝑇sat]=[𝐾].

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 Figure 2.

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 temperature 𝑇am=293K; heat transfer coefficient for free convection 𝛼out=𝐶hot. On the wall inside surface the gas mixture temperature 𝑇mix(𝑡) 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 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 𝐶steel=460J/(kg⋅K) and 𝐶concrete=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 Figure 2) 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. Figure 5 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.

It follows from Figure 6 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 Δ𝑇≪1𝐾. 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: 𝑇=𝑇am.

The simulation results shown in Figure 5 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 𝑡 equals𝑄(𝑡)=𝑥2𝑥1𝜌𝐶𝑇(𝑥,𝑡)−𝑇𝑥,𝑡0𝑑𝑥,(51)
where 𝑇(𝑥,𝑡) is the function of through-the-thickness temperature distribution in the region [𝑥1;𝑥2] with uniform physical properties 𝜌 and 𝐶; 𝑡0 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 times 𝑡∈(0;10000s] for the problem-specific fixed values of parameters (𝑡0=0; 𝑥1=0, 𝑥2=0.006m for steel liner; 𝑥1=0.006m, 𝑥2=0.3m for concrete wall; 𝜌 and 𝐶 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 Figure 7), 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 (𝑡=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 of 𝑙𝑒≤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 Figures 5 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 Figure 8.

Figure 8: Calculated results of the temperature field (K) (a) and isosurface of 40-percent molar fraction of steam (b) in containment rooms.

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 Figure 9.

Figure 9: Integral operation (time) of hydrogen (a) and steam (b).

6. About Models Validation

The examples of 3D CFD analysis results are presented in Figures 4 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 (𝑘−𝜔) 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].

References

M. Babic and I. Kljenak, “CFD spray simulations for nuclear reactor safety applications with lagrangian approach for droplet modelling,” in Proceeding of the International Conference Nuclear Energy for New Europe, p. 13, Portoroz, Slovenia, September 2007.

D. Laerence, “‘Best’ practical guidelines for the use of CFD in nuclear reactor safety applications,” in Proceedings of the CFD Workshop on Test Cases, Databases & BPG for Nuclear Power Plants, p. 18, University of Manchester, July 2008.