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

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∗)

𝜕𝜌𝑌_H2𝜕𝑡+∇⃗⋅𝜌𝑌^H2𝐕⃗=∇⃗⋅𝜇+𝜌𝜈^𝑇Sc_H2∇⃗𝑌_H2−𝜔^H2𝜕𝜌𝑌_O2,(2)𝜕𝑡+∇⃗⋅𝜌𝑌^O2𝐕⃗=∇⃗⋅𝜇+𝜌𝜈^𝑇Sc_O2∇⃗𝑌_O2−𝜉^Omass2𝜕𝜌𝑌_N2𝜔_H2,(3)𝜕𝑡+∇⃗⋅𝜌𝑌^N2𝐕⃗=∇⃗⋅𝜇+𝜌𝜈^𝑇Sc_N2∇⃗𝑌_N2𝜕𝜌𝑌_{atm_H2O},(4)𝜕𝑡+∇⃗⋅𝜌𝑌^atm_H2O=∇⃗⋅𝜇+𝜌𝜈_𝑇𝐕⃗Sc_{atm_H2O}∇⃗𝑌_{atm_H2O},(5)

𝜕𝜌𝑌_H2O𝜕𝑡+∇⃗⋅𝜌𝑌^H2O𝐕⃗=∇⃗⋅𝜇+𝜌𝜈^𝑇Sc_H2O∇⃗𝑌_H2O,(6∗)

𝑌_{PAR_prod}=1−𝑌_H2−𝑌_O2−𝑌_N2−𝑌_{atm_H2O}−𝑌_H2O=1−^𝑁−1_𝑚₌₁𝑌_𝑚,(7)

𝜕𝜌𝐕⃗𝜕𝑡+∇⃗⋅𝜌𝐕⃗𝐕⃗=+𝜌∇−⃗⋅𝜌𝜇_atm+𝐠𝜌⃗−𝜈_𝑇∇⃗𝑃𝜇𝝉−23∇⃗(𝜌𝐾),(8∗)

𝜕(𝜌𝐻)𝜕𝑡=𝜕+𝑃∇⃗⋅𝜌𝐻𝐕⃗+𝜕∇𝑡⃗⋅+𝜕𝜕𝜆𝑄𝑡++𝜆𝜌_𝑇𝐠⃗⋅𝐕⃗+∇⃗𝑇+𝜇+𝜇𝜌𝜈^𝑇𝝉⋅𝐕⃗−23𝜌𝐾_𝑁𝐕⃗_𝑚₌₁∇⃗⋅𝐶_𝑝_𝑚𝑇𝜇+𝜌𝜈^𝑇Sc_𝑚∇⃗𝑌_𝑚,(9∗)

𝜕(𝜌𝐾)𝜕𝑡+∇⃗⋅𝜌𝐾𝐕⃗=∇⃗⋅𝜇+Pr𝜌_𝐾𝜈^𝑇+𝐺−𝛽_∗∇⃗𝐾𝜌𝐾𝜔−𝜈^𝑇Pr_𝑇𝜕(𝜌𝜔)𝐠⃗⋅∇⃗𝜌,(10)𝜕𝑡=∇+⃗⋅∇⃗⋅𝜇+𝜌𝜔𝜌𝐕𝜈_𝑇⃗Pr_𝜔+𝜔∇⃗𝜔𝐾𝛼𝐺−𝛽^∗∗𝜌𝐾𝜔+𝜈^𝑇Pr_𝑇𝑃=𝜌𝑅₀×(𝛼+1)max𝐠⃗⋅∇⃗𝜌;0−𝐠⃗⋅∇⃗𝜌,(11)𝑇^𝑁_𝑚₌₁𝑌_𝑚𝑀_𝑚𝐻=𝐶_𝑝+𝜌_atm𝑔𝑥₃−𝑥_3,0,𝛼_𝑆𝑇+0.5𝐕⃗⋅𝐕⃗;=^𝑁_𝑚₌₁𝜇=𝑓𝑇_𝛼𝑌_𝑚,𝛼_𝑚,𝛼≡𝜇,𝜆,𝐶_𝑝,(12)273.15^3/2273.15+𝐶_𝑆𝑇+𝐶_𝑆𝜆_𝑇𝜇₀,=𝐶^𝑝𝜌𝜈_𝑇Pr_𝑇;Sc_𝑚=𝜇𝜌𝐷_𝑚𝜏,(13)_𝑖𝑗=𝜇𝜕𝑢^𝑖𝜕𝑥_𝑗+𝜕𝑢^𝑗𝜕𝑥_𝑖−32𝛿^𝑖𝑗𝜕𝑢_𝑘𝐺=𝜌𝜈_𝑇𝜕𝑥_𝑘,(14)12𝜕𝜕𝑥𝑢^𝑖_𝑗+𝜕𝑢^𝑗𝜕𝑥_𝑖²−2−23𝜕𝜕𝑥𝑢^𝑘_𝑘²3𝜌𝐾𝜕𝜕𝑥𝑢^𝑘_𝑘𝜈_𝑇,=𝛼¹𝐾𝑆=max𝛼₁𝜔;𝑆𝐹₂,𝑆_𝑖𝑗2𝑆_𝑖𝑗𝑆_𝑖𝑗,=12𝜕𝜕𝑥𝑢^𝑖_𝑗+𝜕𝑢^𝑗𝐹₂𝜕𝑥_𝑖,=tanhmax²2√𝐾𝛽^∗𝜔𝑦;5𝑦0²0𝜔𝑣,(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 𝑢₁, 𝑢₂, 𝑢₃; ∇⃗ is nabla; 𝜌, 𝜌_𝑚, 𝜌_atm are density of mixture, 𝑚th mixture component, and atmosphere at some point in space of containment rooms, respectively; 𝑌_H2=𝑌₁, 𝑌_O2=𝑌₂, 𝑌_N2=𝑌₃, 𝑌_{atm_H2O}=𝑌₄, 𝑌_H2O=𝑌₅, 𝑌_{PAR_prod}=𝑌₆, 𝑌_𝑚=𝜌_𝑚/𝜌 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]; 𝜉^Omass₂ 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]; 𝑅₀ is the universal gas constant; 𝑀_𝑚 is molar mass of the 𝑚th component; 𝑥₃ is the vertical coordinate in the Cartesian frame of reference directed away from the earth center (here, 𝑥₁, 𝑥₂, 𝑥₃ are radius vector coordinates of the point in the Cartesian frame of reference); 𝑥_3,0 is the coordinate of the point 𝑥₃ with a defined value of 𝜌_atm; 𝑓_𝛼(⋯) are known semiempirical functions; 𝐶_𝑆 is the Sutherland constant; 𝜇₀ 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; 𝛼₁=5/9 [7–9]; 𝑦 is the distance to the nearest wall. The functions 𝐺, 𝑆, 𝐹₂ 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𝐕⃗=∇⃗⋅𝜇+𝜌𝜈^𝑇Sc_H2O∇⃗𝑌_H2O+𝜔^H2O;(6)

𝜕(𝜌𝐻)𝜕=𝑡𝜕𝑃+∇⃗⋅𝜌𝐻𝐕⃗𝜕𝑡+𝜕𝜕𝑄𝑡+0.5𝑄^{phase_transition}+×∇𝜔⃗⋅_H2O𝜆−+||𝜔𝜆_H2O𝑇||+𝐶_𝑝_H2O𝑇_sat𝜔_H2O+𝜌𝐠⃗⋅𝐕⃗+∇⃗𝑇+𝜇+𝜇𝜌𝜈^𝑇𝝉⋅𝐕⃗−23𝜌𝐾_𝑁𝐕⃗_𝑚∑₌₁∇⃗⋅𝐶_𝑝_𝑚𝑇𝜇+𝜌𝜈^𝑇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: 𝐽=𝐽₁=𝜌𝐷^steam0.5𝑑_dropln1−𝑌^steam,∞1−𝑌_steam,sat,(16)𝐽=𝐽₂=𝜆^steam0.5𝑑_drop𝐶_𝑝_H2Oln1+𝐶^𝑝_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=𝜋𝑑²_drop, 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𝐶_𝑝_H2Oln1+𝐶^𝑝_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𝜋𝑑³_drop6,(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/3_drop⋅𝑑𝑚^drop=−𝜔_H2O𝑑𝑡3𝑛𝜋𝜌6_liquid^1/3𝑚^−2/3_drop.(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𝑑𝑚_dropor𝑑𝑑𝑡𝑑𝑡𝐻^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}=−𝜔_H2O𝐶_liquid𝑄_{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.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=𝑛𝜋𝑑²_drop𝐽=2𝑛𝜋𝑑_drop𝜌𝐷_steamln1−𝑌^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,satsteam,sat}𝑇^sat,𝑌_𝑚,𝑚=1,6𝑇^drop,𝑌_𝑚,𝑚=1,6,(27) the above equation for 𝜔_H2O will be written as 𝜔_H2O=2𝑛𝜋𝑑_drop𝜌𝐷_steam×ln⎡⎢⎢⎣1−𝑌^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=𝜋𝑑²_drop/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⁻⁷𝑢𝑑𝜈^drop+0.466+45.822𝑢𝑑𝜈^drop⁻¹−28.262𝑢𝑑^drop𝜈⁻²+11.995𝑢𝑑^drop𝜈⁻³,(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𝐠⃗−𝜋𝐶^𝑥𝜌𝑑²_drop𝐕_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, 𝑉₀. It contains (𝑛𝑉₀) droplets. Forces acting from them are given by the following summand:−_𝑆drop𝝉⃗_𝑛𝑑𝑆𝑛𝑉₀=𝑛𝑉₀𝜋𝐶_𝑥𝜌𝑑²_drop𝐕_drop||𝐕_drop||8.(34) Dividing the summand by 𝑉₀ and letting 𝑉₀ 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:𝑛𝜋𝐶_𝑥𝜌𝑑²_drop𝐕_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⋅𝜌𝐕⃗𝐕⃗−2𝐠⃗−∇⃗𝑃+∇⃗⋅𝜇+𝜇𝜌𝜈^𝑇𝝉3∇⃗(𝜌𝐾)+𝑛𝜋𝐶^𝑥𝜌𝑑²_drop𝐕_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 𝑉₀. 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=𝑆⁻¹_drop∬_𝑆drop𝝉⃗^drop_𝑛𝑑𝑆. Then,_𝑆drop𝝉⋅𝐧_mix⋅𝐕_drop𝑑𝑆=≈−−𝝉⃗_𝑆drop𝝉_{drop𝑛,aver}⃗^drop_𝑛⋅𝐕_drop𝑑𝑆⋅𝐕_drop𝑆_drop==−𝜋𝐕C_drop𝑥⋅_𝑆drop𝝉⃗^drop_𝑛𝑑𝑆𝜌𝑑²_drop||𝐕_drop||³8.(36) The small volume 𝑉₀ contains (𝑛𝑉₀) droplets; consequently, the summand characterizing the forces exerted by the droplets in the energy equation for the volume 𝑉₀ will be given by the expression: 𝑛𝑉₀𝜋𝐶_𝑥𝜌𝑑²_drop|𝐕_drop|³/8. We divide the summand by 𝑉₀ and let it go to zero: 𝑛𝜋𝐶_𝑥𝜌𝑑²_drop|𝐕_drop|³/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+𝜌𝐠⃗⋅𝐕⃗+∇⃗𝑇+𝜇+𝜇𝜌𝜈^𝑇𝝉⋅𝐕⃗−23𝜌𝐾_𝑁𝐕⃗_𝑚₌₁∇⃗⋅𝐶_𝑝_𝑚𝑇𝜇+𝜌𝜈^𝑇Sc_𝑚∇⃗𝑌_𝑚+𝑛𝜋C_𝑥𝜌𝑑²_drop||𝐕_drop||³8.(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)

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]: 𝑗_𝑏=𝑞^wall𝑄_{phase_transition}𝑞_wall,(𝑥)=𝛼(𝑥)𝑇_sat−𝑇_wall,𝑞_wall𝛼𝐻(𝑥)==𝛼𝜆𝐻_3liquid𝑇_sat−𝑇_wall,𝜌²_liquid𝑔𝑄_{phase_transition}4𝜇_liquid𝑇_sat−𝑇_wall𝑥^0.25,𝛼𝐻=43𝜆^3liquid𝜌²_liquid𝑔𝑄_{phase_transition}4𝜇_liquid𝑇_sat−𝑇_wall𝐻^0.25𝛿(𝑥)=4𝜆_liquid,𝜇_liquid𝑇_sat−𝑇_wall𝑥𝑄_{phase_transition}𝜌²_liquid𝑔^0.25,𝑤(𝑥)=𝜌^liquid𝑔3𝜇_liquid𝛿²(𝑥),(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<Re_crit=400 [19], where Re is the (film) Reynolds numbers: Re=𝜌_liquid𝑤𝛿/𝜇_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, |𝐕⃗|≤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 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=𝑞_wall1+𝑄𝑞^steam_{phase_transition}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 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 pressure 𝑝₀ of the steam-air mixture is the sum of partial pressures 𝑝_steam and 𝑝_air of steam and air, respectively, that is, 𝑝₀=𝑝_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<𝑝₀, the temperature 𝑇_sat is always lower than the saturation point at pressure 𝑝₀.
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,𝑝,𝑇3⁴⎷𝜆³_liquid𝜌²_liquid𝑔𝑄_{phase_transition}×14+𝜇𝐶_liquidsteam𝑇_sat(𝑝)−𝑇_wall𝐻𝑇−𝑇_sat(𝑝)𝑄_{phase_transition}𝑇^sat(𝑝)−𝑇_wall𝑇−𝑇_wallif𝑇_sat,0,if𝑇_sat(𝑝)>𝑇_wall;(𝑝)≤𝑇_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.16i𝛼f0_hot.1𝐻6,𝑟^H₂^O𝑝,𝑇,𝐶_hot𝛼_hot𝐻,𝑟^H₂^O𝑝,𝑇>𝐶_hot;i,f0.16𝛼_hot𝐻,𝑟^H₂^O𝑝,𝑇≤𝐶_hot,(43) where 𝐶_hot≈5.9W/(m²K). 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}⁻¹.(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.32exp18.76−4081.18𝑇_sat−36.38(49) or 𝑇_sat(𝑝)=36.38−4081.18ln𝑝133.32−18.76⁻¹,(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 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/m³ and 𝜌_concrete=2400kg/m³specific 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.

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 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𝜌𝐶𝑇(𝑥,𝑡)−𝑇𝑥,𝑡₀𝑑𝑥,(51) where 𝑇(𝑥,𝑡) is the function of through-the-thickness temperature distribution in the region [𝑥₁;𝑥₂] with uniform physical properties 𝜌 and 𝐶; 𝑡₀ 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, 𝑥₂=0.006m for steel liner; 𝑥₁=0.006m, 𝑥₂=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.

(a)

(b)

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.

(a)

(b)

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].