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

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.


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][2][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][5][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][5][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 steelliner containment walls accounting for steam condensation on them.The spray system is installed under the dome of the central reactor building.According to its purpose, it belongs to confinement safety systems, and its basic functions include pressure and temperature control in the containment during severe accidents.
Combustion risk of such clouds in this case is assessed using the known Shapiro diagram.In the first approximation, the aqueous solution of boric acid distributed by the spray system can be considered similar to water in its physical properties.
Within our model we assume that the gas mixture in containment rooms has six components: hydrogen produced during severe accidents, (m = 1); air as a mixture of its three basic components-oxygen, (m = 2), nitrogen, (m = 3), and atmospheric water vapor, (m = 4); steam injected into the containment during severe accidents, (m = 5); additional steam produced by hydrogen and oxygen absorption by PAR [4][5][6] (referred to as absorption product) (m = 6).

Selection of Underlying Model Equations
As a basis for the development of the fluid dynamics model of processes inside the containment, we used an adapted set of Reynolds equations completed by the (k − ω)-turbulence model [7][8][9].
The turbulence model was chosen by the authors based on the results of computational experiments with steamhydrogen-air flows in containment rooms without considering the operation of the spray system.These experiments demonstrated considerable advantages of this model over the zero-order turbulence models and various modifications of the (k − ε)-turbulence model [10,11].Similar conclusions have been made in [12,13].Application of the zero-order turbulence models in most of computational experiments lead to clearly physically invalid numerical results for the space-time field of gas flow velocities in the simulated rooms.
Thus, the basis for the fluid dynamics model of processes inside the containment with neglect of the spray system operation can be represented in the following way (see, e.g., [7][8][9]14]): where ) is the substantial derivative of the scalar function (writing a substantial derivative of a vector function means substantial differentiation of vector function components); t is time; V is velocity with components u 1 , u 2 , u 3 ; ∇ is nabla; ρ, ρ m , ρ atm are density of mixture, mth mixture component, and atmosphere at some point in space of containment rooms, respectively; ρ is the relative mass fraction of the mth component; μ is the dynamic viscosity; μ T = ρν T is the turbulent viscosity; ν is the kinematic viscosity; Sc m is the Schmidt number corresponding to the mth component; ω H2 is the mass rate of the generalized irreversible single-stage reaction of hydrogen absorption in PAR [4,5]; ξ mass O2 is the stoichiometric mass coefficient of oxygen; N is the number of components in the gas mixture (in our case, N = 6 (see above)); the subscript T stands for "turbulent," and the subscript atm, for the "atmosphere in containment rooms"; g is gravity acceleration (g is the gravity acceleration modulus); P is piezometric pressure of mixture; τ is tensor of viscous stresses with components τ i j ; K is turbulent kinetic energy; H = h + 0.5 V • V is total specific (per unit mass) enthalpy (h = C p T is enthalpy for perfect gas, C p is specific (per unit mass) heat capacity at constant pressure; T is mixture temperature); ∂Q/∂t is the specified rate of heat release from outer sources per unit volume (heat generation rate per unit volume); λ is the thermal conductivity; λ T is the turbulent thermal conductivity; Pr is the Prandtl number (the subscript K indicates that the value of the Prandtl number is defined specifically for turbulent energy equation (10), and ω indicates that it is defined for turbulence dissipation (11)); Pr K = 2.0 [7][8][9][10][11]; β * = 0.09 [7][8][9]; ω are turbulence frequencies; Pr ω = 2 [7][8][9]; α = 5/9 [7][8][9]; β * * = 0.075 [7][8][9]; R 0 is the universal gas constant; M m is molar mass of the mth component; x 3 is the vertical coordinate in the Cartesian frame of reference directed away from the earth center (here, x 1 , x 2 , x 3 are radius vector coordinates of the point in the Cartesian frame of reference); x 3,0 is the coordinate of the point x 3 with a defined value of ρ atm ; f α (• • • ) are known semiempirical functions; C S is the Sutherland constant; μ 0 is the dynamic viscosity under normal conditions; D m is the coefficient of binary diffusion of the mth component in the rest of the mixture; δ i j is the Kronecker delta; α 1 = 5/9 [7][8][9]; y is the distance to the nearest wall.The functions G, S, F 2 in (15) are auxiliary.Note that the model uses Favreaveraged velocity components and thermal variables (H, h, T) [15] (i.e., mixture density is used as a weighting function) and classical Reynolds-averaged density and pressure [11].Asterisks in the right upper corner in (1)- (15) and in the following indicate the equations that will be modified as the model will be constructed.
It should be noted that in the model ( 1)-( 15) steam was conventionally divided into 2 components, such as atmospheric steam and injected (during accident) steam.That was done for the implementation of comprehensive analysis of distribution parameters of injected (during accident) mixture in containment rooms.

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: 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; Q phase transition is specified latent mass heat of steam condensation into water droplets at ω H2O < 0; T 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 Q 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.5Q phase transition [ω H2O − |ω H2O |].The summand in (9 * * ) given by [(C p ) H2O T 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 J from a spherical droplet: where ρ is gas mixture density; D steam is the coefficient of steam diffusion in the mixture; d drop is droplet diameter; Y steam,∞ and T ∞ are the mass fraction and temperature of steam far from the surface, respectively (where effects produced by the evaporating steam can be ignored); Y 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 S drop = πd 2 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: where dm steam drop /dt is the rate of steam production from one droplet; m drop is droplet mass; dm drop /dt is the rate of droplet mass variation.Let n be a specific (per unit volume of matter) number of spherical water droplets.In this case, the mass of steam produced from all unit-volume droplets can be calculated using the formula (see (17)) Here, as with ( 1)-( 15), the temperature T is equal to the gas mixture temperature (far from the droplet surface).
As it evaporates, the droplet decreases in its diameter.Let us estimate the rate of droplet diameter variation.The droplet mass is where ρ liquid is density of the droplet liquid (water); V drop is droplet volume.Hence, d drop = [6m drop /(πρ liquid )] 1/3 , and 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: T sat ≈ T 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), h liquid = ρ liquid C liquid T liquid , where h liquid is specific (per unit volume) enthalpy of the liquid; ρ liquid = const is liquid density (the liquid is incompressible); C liquid = const is specific (per unit mass) heat capacity of the liquid; T liquid is liquid temperature.Considering this EOS, droplet enthalpy can be written as where H drop and T drop are average enthalpy and temperature of the droplet, respectively.Hence, 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 T sat : By equating the right-hand sides of the second equation in ( 22) and ( 23), we have 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 for the following boundary conditions: , where r is distance from the droplet center to an arbitrary point.
During droplet evaporation, Y steam,∞ ≤ Y steam,sat , which produces a flow of steam from the droplet surface to the gas mixture.If the situation is opposite, Y steam,∞ > Y steam,sat (i.e., the steam mass fraction in the gas mixture is higher than on the droplet surface equal to Y steam,sat ), a reverse flow of steam-from the gas mixture to the droplet surfaceshould 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, If we first write the mass fraction of steam on the saturation line Y steam,sat as a function of temperature (for known mass fractions of gas mixture components far from the droplet), the above equation for ω H2O will be written as In order to estimate Y steam,sat for the complex (T sat , {Y m , m = 1, 6}), one can employ the following algorithm: based on the known value of T 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 V 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)): where p is static pressure of the gas mixture.This equation is integrated over the droplet volume.We assume here that liquid density is constant, and pressure variation in the droplet volume can be ignored: Science and Technology of Nuclear Installations Droplet drag is estimated here using a function that approximates experimental data [14]: where C x is the drag coefficient of the droplet; ρ is density of the gas mixture around the droplet; S MS = πd 2 drop /4 is midsection area.Here, just as in (30), τ n is the projection of τ to the outward normal to the droplet surface.Tabulated reference parameters of the function C x = C x (u) (where u is the velocity of droplet motion relative to the gas environment) have been approximated by Kobyakov in [14]: where ν is the kinematic viscosity of the gas mixture.In the first approximation, we assume that Thus, the equation of droplet motion can be written in the following way: 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, V 0 .It contains (nV 0 ) droplets.Forces acting from them are given by the following summand: Dividing the summand by V 0 and letting V 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: 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 * ): Let us now proceed to transformations of energy equation (9 * ).As before, let us consider a small volume of the gas mixture V 0 .The power of forces experienced by the gas mixture as a result of droplet motion equals where n mix is unit vector normal to the droplet surface external to the gas mixture; n drop is unit vector normal to the droplet surface external to the droplet.Let τ drop n = τ • n drop .We introduce the notion of droplet surface average viscous stress: The small volume V 0 contains (nV 0 ) droplets; consequently, the summand characterizing the forces exerted by the droplets in the energy equation for the volume V 0 will be given by the expression: We divide the summand by V 0 and let it go to zero: To write the final form of the energy equation, one needs to add the resulting summand to the right-hand side of (9 * ): The velocity of droplet motion is V drop .Consequently, it can be assumed that the variable n (specific (per unit volume of gas mixture) number of droplets) also moves at a velocity of V drop .Based on these considerations, one can write an equation for n: Thus, the mathematical model of spray system operation will combine (1)-( 15), ( 20), ( 21), ( 24), ( 28), ( 32), ( 33), (37)).

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][20][21][22]: , where j b is the mass flux of condensing steam; q wall is the heat flux during steam condensation; x is the distance from an arbitrary point on the wall to the upper film surface (vertically); H is the film height; q wall (x) and α(x) are local (for the coordinate x) values of the heat flux and heat transfer coefficient; q wall ( H) and α( H) are film average (for the film height H) values of the heat flux density and heat transfer coefficient; λ liquid , ρ liquid , and μ liquid are the thermal conductivity, density, and dynamic viscosity of condensate (as a rule, these are considered for the temperature 0.5(T sat + T wall )); T wall is the wall temperature; δ(x) is the local thickness of the condensate film; w(x) is the velocity of condensate averaged through the film thickness.The wall temperature, saturation point, and liquid properties are assumed constant.Laminar flow of the film is observed subject to the condition Re < Re crit = 400 [19], where Re is the (film) Reynolds numbers: Re = ρ liquid wδ/μ liquid .
Effects of Steam Flow Velocity.Composite formula (38) has limited application, because it has been developed for fixed dry saturated steam.However, in accordance with [16], this formula remains valid for small velocities of steam flow, | V| ≤ 10 m/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 q steam = C steam (T steam − T sat ) and replace evaporation heat Q phase transition in the formula by heat of superheating Q phase transition = Q phase transition + q steam , where T steam is the temperature of superheated steam, C steam is heat capacity of superheated steam [16].Temperature difference is again assumed to be equal to ΔT = T sat − T 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 Q phase transition > Q 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)) q wall = α(T sat − T wall ), j b = q wall /Q phase transition .
If steam is superheated, additional heat flux, q steam add = q steam j b = q steam q wall /Q phase transition , will be withdrawn from the gas mixture.Thus, the total heat flux toward the wall will be q wall,Σ = q wall + q steam add = q wall 1 + q steam Q phase transition or where α Σ is the wall-to-superheated-steam heat transfer coefficient; note that here T steam = T. Hence, 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 p 0 of the steam-air mixture is the sum of partial pressures p steam and p air of steam and air, respectively, that is, p 0 = p steam + p air .As a result of steam condensation, the pressure p steam near the wall is lower than in the rest of the volume.Therefore, the pressure p steam decreases continuously as steam approaches the wall, and the closer it comes to the wall, the faster it decreases, whereas the pressure p air increases [16].This produces high air mass fraction near the wall, which is stands for "air"; subscript "Π" stands for "steam") [16].
identified as a layer that can be penetrated by steam molecules only by way of diffusion.
As a result of these processes, the temperature difference ΔT = T sat − T wall also decreases, because the temperature of the mixture is always equal to the temperature of steam saturation at partial pressure p steam .Since p steam < p 0 , the temperature T sat is always lower than the saturation point at pressure p 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 r H2O be the volume fraction of steam in the mixture.Then, the partial pressure of steam p H2O can be estimated as p H2O = r H2O p.The saturation point of steam can be calculated in this case as T sat = T sat (p H2O ) = T sat (r H2O p), where T sat = T sat (p) is the temperature as a function of pressure on the saturation line of water.
For numerical simulations of steam condensation on the inner walls of a reactor building, we introduce a parameter α hot to characterize the average heat transfer coefficient of single-component superheated resting steam on a vertical flat wall (see (38)): The condition in ( 42) is introduced to constrain the range of application for the steam condensation formula.In particular, for T sat < T 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 where C hot ≈ 5.9W/(m 2 K).The value of the heat transfer coefficient due to free convection (without condensation processes) equals α = C 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 q wall = α(T − T 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 T to T sat ), the mass flux of condensation can be calculated using the formula 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: where ρ wall , C wall , T 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: Here and below, n is unit outward normal to the surface of the region of interest.On the gas mixture side of the same interface one should use the Neumann boundary conditions for steam "withdrawal": By analogy with (47), the boundary condition for the energy equation will be to characterize the heat flux generated by gas-to-liquid "transition" of water (for the mass flux j b ) at T sat .
The saturation point of steam can be calculated using the equation similar to the Antoine equation [23]: where dimensionality in the formula has the form of

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.
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 T am = 293 K; heat transfer coefficient for free convection α out = C hot .On the wall inside surface the gas mixture temperature T mix (t) varied by the law, depicted in Figure 3.For the heat transfer coefficient was taken equality α in = α out .Such a value of α in was chosen based on the results of preliminary CFD analysis of steam-hydrogen mixture flow in reactor building rooms (Figure 4).Physical properties of the materials were specified in accordance with [27]: steel and concrete densities ρ steel = 7850 kg/m 3 and ρ concrete = 2400 kg/m 3 specific heat capacities of steel and concrete C steel = 460 J/(kg • K) and C concrete = 780 J/(kg • K); thermal conductivities of steel and concrete λ steel = 45 W/(m • K) and λ concrete = 1.23 W/(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.
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 ÷ 300 mm.In the wall's concrete layers situated far from the heated surface, temperature variation during the accident stays within ΔT 1 K. Therefore, for the numerical analysis of severe accidents in the containment  wall/environment heat exchange it would suffice to supplement the computational model with solid domain consists of steel liner and 250 mm concrete layer.On the outside surface of the concrete layer one should assign the Dirichlet boundary condition: T = T 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 t equals where T(x, t) is the function of through-the-thickness temperature distribution in the region [x 1 ; x 2 ] with uniform physical properties ρ and C; t 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 t ∈  (0; 10000 s] for the problem-specific fixed values of parameters (t 0 = 0; x 1 = 0, x 2 = 0.006 m for steel liner; x 1 = 0.006 m, x 2 = 0.3 m for concrete wall; ρ and C corresponding to the above values for steel and concrete) gives time profiles of heat accumulated in the wall's structural parts.These profiles are shown in Figure 7.
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 (t = 10000 s) 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 l e ≤ 0.006 m 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.
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.

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 (k − ω) turbulence model in (10,11) with regard to steam flow was done in comparison of the results of numerical analysis with data from [29][30][31][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][4][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 CFDanalysis of the emergency emissions of steam-hydrogen-air mixture inside containment.

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

Figure 1 :
Figure1: 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].

Figure 2 :Figure 3 :
Figure 2: Sketch of a model heat transfer problem to estimate the parameters of heat exchange through the inner containment wall.

1 )Figure 4 :
Figure 4: Example of a snapshot of the field of velocities (m/s) in containment rooms.

Figure 5 :Figure 6 :
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 7 :
Figure 7: Heat accumulated in the wall's structural parts during the severe accident.

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.