The article gives a general description of the SIMCO calculation code designed to simulate thermohydraulic and physico-chemical processes in containments of nuclear power facilities. The authors present a calculation technique based on a physico-mathematical model in lumped parameters. As a numerical solution method, the modified semi-implicit SIMPLER procedure is used. The code was examined using analytical and qualitative tests. A comparison of the numerical and analytical solutions showed good agreement. The code was verified using the experimental data obtained at the NUPEC installation (Japan). Based on the results of testing and verification, it was concluded that, in general, physico-mathematical code models adequately describe the processes of heat/mass transfer in the containment. Therefore, this SIMCO code version can be used to analyze the totality of thermophysical and physico-chemical processes in nuclear power facilities with containments, including the transfer of hydrogen/steam/air mixtures.

SIMCO code, physico-mathematical model, lumped parameters, analytical test, verification, hydrogen distribution

Modernization of existing and development of new nuclear power plants meeting the increased requirements for reliability, safety and economical efficiency make it necessary to improve the calculation methods for reactors, heat exchange equipment as well as safety and accident localization systems. The containment is the fourth, final safety barrier to the spread of radioactive products to the environment. It is designed in accordance with regulatory documents, and heat/mass transfer processes in the containment during loss-of-coolant accidents (LOCA) are complex spatial in nature: they are characterized by numerous thermal, physical and chemical phenomena that determine the safety system effectiveness and reliability. For this reason, research processes in the containment become an integral part of design works. For new generation NPPs, special attention is paid to hydrogen fire/explosion safety (IAEA-TECDOC-1196 2001, Kirillov et al. 2017, Shapiro and Moffette 1957, Soldatov and Golodnova 2009). The urgency of the problem of efficient and reliable hydrogen explosion safety increased after the accident at the Fukushima Daiichi NPP (Japan) on March 11, 2011.

To simulate containment of NPPs, it is possible to use the SIMCO code originally developed as the KUPOL-M code real-time version (Kirillov et al. 2017). Simulating dynamic processes in real time (or as close as possible to this) is particularly topical in developing analytical and full-scale simulators of existing and projected nuclear power facilities. The presence of simulators at nuclear power facilities provides the most efficient training of NPP operating personnel. The SIMCO code is designed to calculate the thermophysical and physico-chemical parameters of the containment environment for nuclear power plants or other nuclear power facilities with containments under various accident scenarios.

Simulating is performed using an arbitrary topology of control volumes inside the containment area.

The following basic values are calculated:

• time changes of temperature and gas pressure in the containment;
• non-stationary temperature distribution in the walls and equipment;
• time-density dependences of nitrogen, oxygen, steam, hydrogen, helium, carbon dioxide, carbon monoxide, an arbitrary inert component, and fine moisture in the containment;
• steam condensation intensity and condensate temperature in each compartment;
• temperature of water drained to the lower compartment trays;
• gas mixture flowrate in penetrations between the compartments.

The calculation takes into account the effects of unsteady gas mixture heat/mass transfer, steam volume/surface condensation, gas mixture natural convection as well as the operation of the sprinkler system and passive catalytic hydrogen recombiners, pumps and ventilation systems.

Applicability constraints:

1. The use of the SIMCO code is limited to the region of atmospheric thermodynamic parameters in the containment in case of LOCAs at NPPs.
2. The code requires the following boundary conditions:
   • flowrate and enthalpy (temperature) of the steam-water mixture entering under the containment from the leak;
   • flowrates and temperatures of steam-gas components entering the compartments;
   • exhaust ventilation flowrate;
   • flowrate and enthalpy of water supplied to the sprinklers.
3. Acceptable values of the gas-mixture parameters:
   • pressure from 0.07 to 1.5 MPa;
   • temperature from 275 to 1200 K.

The mathematical model of heat/mass transfer in the containment (Zaitsev 2005, Levchenko 2007) is a system of ordinary differential equations for conservation of momentum, energy and gas mixture components recorded for each selected control volume. This approach (i.e., modeling in lumped parameters) has become a frequent practice (Gauntt et al. 2001, Murata 1997, Povilaitis et al. 2015). In view of the significant effect of compressibility on the gas flow characteristics, instead of the momentum-conservation equation, dependences obtained for the adiabatic gas flow from pressure vessels are used. The universal gas law has been adopted as the equation of state. The model of heat exchange with walls and various equipment includes a one-dimensional heat conduction equation for flat, cylindrical, and spherical walls. Empirical dependences are simulated for the total heat transfer coefficient (to calculate the surface steam condensation) and for convective heat transfer coefficients. The process of volumetric steam condensation includes modeling droplet-size spectra.

The numerical solution to the obtained system of equations is carried out on the basis of the modified semi-implicit procedure SIMPLER (Patankar 1984, Launder 1976, Samarsky and Nikolaev 1978, Bubenchikov and Kharlamov 2001). The scalar values (pressure, temperature and density) are determined in the control volumes, and the steam-gas-droplet mixture velosi are defined at the boundaries of the control volumes. At the first stage, the pressure in the control volumes and the mixture flowrates in the bonds between them are determined. Next, the equations of energy transfer and mixture components are solved. At the last stage of the time step calculation, local problems of volume and surface steam condensation are solved and temperature fields in the walls and equipment are determined.

The finite-difference system of algebraic equations, which generally has a sparse matrix of coefficients, is directly solved by the Gauss method or the lower relaxation method. The differential equation for the steam-gas mixture flowrate in the bonds is solved by the Runge-Kutta method. The heat equation for the walls, approximated by an implicit difference scheme, is solved by the sweep method.

Analytical “gas-dynamic spring” test. A similar problem for the linearized flow was considered in (Landau and Lifshits 1988). Also, a similar test was used to test the KUPOL-M code (KUPOL-M code. Ver. 1.1 2002).

Let us consider the analytical solution of the linear problem of gas compression in an open vessel under external constant pressure. In order to linearize the problem, we neglect the gas temperature change due to the external influence, then the equation of motion can be written as

(1)

where u is the gas velocity in the bond; ξ − is the linear hydraulic resistance coefficient, m/s; r₀, P₀ is the gas density and pressure in the external environment (kg/m³ and Pa, respectively); DX is the characteristic bond length, m; Р is the pressure in an open vessel, Pa.

The gas density equation in an open vessel is as follows

(2)

where S is the bond cross section; V is the vessel volume.

Using the universal gas law, we replace the pressure in the vessel in equation (1) with the density and, after time differentiation, obtain

ü + au + bu = 0, (3)

where a = ξ/(2DX); b = RT₀S/(M_GVDX); R is the universal gas constant; T₀ is the gas temperature in the vessel; M_G is the gas molecular weight.

Assuming that the gas in the vessel is identical to the external gas and considering the initial conditions

(4)

we obtain the following solution to the equation

u = 2DP₀λ^–1 exp(aτ) sin(λτ/2), (5)

where λ = (4b – a²)^1/2; DP₀ is the initial pressure drop between the boxes.

The numerical simulation of this test was carried out using a two-box nodalization scheme with a linearized model for the equation of motion in the isothermal approximation. The external environment is simulated by the larger box, the volume of which is several orders of magnitude greater than the that of the open vessel. With equal initial gas densities in the boxes, the required pressure drop is created due to the weight of the gas column in the larger box.

The equation of gas density in an open vessel (2) is valid only for positive rate values (gas enters the open vessel from the environment), since the environmental gas density in the right hand side of this equation is constant. Therefore, the numerical solution is compared with the analytical solution in the initial period of time. The geometric characteristics of the nodalization scheme are given in Tab. 1.

Table 1.Download as CSV 

The input data for the “gas-dynamic spring” test

Upper box volume, m3	106
Lower box volume, m3	1.0
Upper box height, m	10.0
Bond flow cross section, m2	0.1
Characteristic bond length, m	1.0
Hydraulic resistance coefficient, m/s	0.5
Initial gas density, kg/m3	1.0
Temperature, K	333
Gas molecular mass, kg/mol	0.028

The dependence of the absolute solution error is shown in Fig. 1 (maximum relative error is ~0.8%). Note that the error depends almost linearly on the time integration step of the motion equation.

Figure 1. 

Difference between the analytical and numerical solutions

The numerical solution after the time point τ > 0.03 s is oscillatory. Damping of oscillations occurs as a result of the mechanical power dissipation due to hydraulic resistance. In real objects described by the square resistance law (as well as an increase in temperature due to the work of external compression forces), oscillations damp out much faster.

Qualitative “natural gas circulation” test. This test illustrates the probability of a gas mixture “parasitic” natural circulation in a closed-loop system if the numerical model does not take into account the gas compressibility due to its own weight throughout the height of the calculated box. This effect is particularly important in modeling concentration/temperature stratifications of the light component (hydrogen) (Gupta et al. 2015).

Let us consider the circulation loop (its nodalization scheme is shown in Fig. 2).

Figure 2. 

Nodalization scheme of the closed “natural circulation” loop

The loop is thermally insulated, there are no heat and mass sources, the gas is homogeneous.

With the height invariable dencity gas inside the box, the pressure throughout its height varies linearly. If the pressure drop across the top bond is zero, then the gas density in both branches at the top elevation is the same and equal to r₀. The pressure in the left branch at the “0 m” elevation is equal to Р₀ + 40r₀g. The pressure in the right branch at the “0 m” elevation will be greater, since the branch height-average gas density due to its compressibility in calculated Boxes 4, 3 and 2 is greater than r₀, i.e., this nodalization scheme (if the Boltzmann density distribution throughout the box height is ignored) leads to the appearance of a non-zero pressure drop, causing the gas flow.

Figures 3 and 4 show the results of numerical calculations of the considered closed loop with parameters close to the actual values of modern NPP containments (Table 2); the calculations were based on two mathematical models (one of height invariable dencity gas inside the box and the other of the SIMCO code, taking into account the gas compressibility (Povilaitis et al. 2015)).

Table 2.Download as CSV 

Parameters of the closed “natural circulation” loop

Box	1	2	3	4	5
Box floor elevation, m	0.0	0.0	10.0	20.0	30.0
Box ceiling elevation, m	40.0	10.0	20.0	30.0	40.0
Box initial gas density, kg/m3	1.14	1.14	1.14	1.14	1.14
Bonds
Bond cross section, m2			0.1
Bond hydraulic resistance coefficient			0.5

Figure 3. 

Occurrence of a stationary “parasitic” natural circulation in the loop (the model ignores the gas compressibility inside the box)

Figure 4. 

Absence of a stationary circulation in the loop (SIMCO code model)

The gas “parasitic” circulation rate in the loop without sources approximates to 0.7 m/s. The SIMCO code model, which takes into account the steam-gas mixture compressibility in the calculated box volume, after the transient stage associated with the initial conditions, shows the absence of a “parasitic” circulation in the system of connected boxes. The transient period duration is determined by the hydraulic characteristics, box geometry and gas mixture composition. Thus, the incorrectness of the physical-mathematical natural gas convection model can lead to significant errors in calculations of circulation loops in the containment, which is especially important when stratified flows are calculated in the hydrogen safety analysis.

The code was verified using the experimental data obtained at the NUPEC installation (Japan). The NUPEC integrated installation is a model of a PWR containment (linear scale 1: 4) with a volume of 1300 m³ (height = 17.4 m; internal diameter = 10.8 m). The inner space of the containment is divided into 25 compartments simulating the interior of the prototype containment.

Table 3 shows the main geometric characteristics of the containment rooms model.

Table 3.Download as CSV 

Rooms characteristics

Cross section, m2	Volume, m3	Floor elevation, m	Height, m
4.54	5.95	2.097	3.316
6.76	14.95	3.200	2.213
13.06	28.91	3.200	2.213
14.84	32.84	3.200	2.213
14.84	32.84	3.200	2.213
13.06	28.91	3.200	2.213
4.88	10.80	3.200	2.213
5.72	12.66	3.200	2.213
1.46	3.74	3.200	2.563
5.72	12.66	3.200	2.213
4.88	10.80	3.200	2.213
28.10	53.05	5.425	1.888
28.10	53.05	5.425	1.888
4.88	9.22	5.425	1.888
5.72	10.80	5.425	1.888
1.46	2.14	5.775	1.463
5.72	10.80	5.425	1.888
4.88	9.22	5.425	1.888
12.70	24.14	5.425	1.900
1.465	2.01	7.325	1.371
1.47	2.016	7.325	1.371
1.46	5.44	7.250.	3.725.
1.47	2.016	7.325	1.371
1.465	2.01	7.325	1.371
91.8	931.3	7.325	12.1175.
Total	1312.272

A series of integral experiments were conducted on this experimental installation. In contrast to the experiments on the HDR installation (Valencia et al. 1992), the NUPEC installation simulated the sprinkler system operation in the containment. The SIMCO code was verified using the M-7-1 experiment (hydrogen mixing and distribution in the containment, the sprinkler system functioning). This experiment was used as the international standard problem ISP-35. For safety purposes, helium was used as a hydrogen simulator. The purpose of the experiment was to study the helium distribution in the containment during the sprinkler system operation. The main factors affecting the distribution of the light component (helium) were natural convection and the sprinkler system operation.

The necessary initial conditions in the compartments of the containment model (1.4 bar, 70 °C) were provided by preliminary heating due to steam supply for 3.5 hours. During the experiment, the steam-helium mixture was fed with a linear change in the flow of the components. The steam flow decreased from 0.08 to 0.03 kg/s during the injection period of 30 minutes. The helium flow during the same period increased from zero to a maximum (0.03 kg/s) and again decreased to zero. The supplied sprinkler water flow was kept constant for 30 minutes being equal to 19.4 kg/s.

According to (Test M-7-1 1994), the following parameters are available for comparison:

1. pressure in the experimental installation;
2. two points on the steam-gas mixture temperature;
3. two points on the light component (helium) local concentrations.

Although the initial conditions in the containment model established as a result of preheating contain a number of uncertainties, it is possible to verify the code using available experimental data.

A compartment with a steam/helium leak was simulated by two calculated volumes.

Figures 5–7 show the results of a comparison of the calculated and experimental dependences of the thermodynamic atmospheric parameters in the compartments of the containment model (pressure, temperature and concentration of helium as a hydrogen imitator).

Figure 5. 

Pressure dynamics in the M-7-1 experiment

Figure 6. 

Helium concentration dynamics in the compartment with steam-helium mixture leak in the M-7-1 experiment

Figure 7. 

Helium concentration dynamics in the dome space in the M-7-1 experiment

The maximum discrepancies between the calculated and experimental data for the gas mixture pressure and temperature are observed at the sprinkler system shutdown stage. Similar results for this period of the experiment were obtained when calculations were made by other containment codes (Test M-7-1 1994) within the international standard problem ISP-35. The true scenario of the experiment may have been somewhat different from that described in (Test M-7-1 1994) at the the sprinkler system shutdown stage. The calculated data on the light component (helium) concentration are in good agreement with the experimental data. The maximum error does not exceed 15% in the compartment with steam-helium mixture leak. It should be considered that this is a local characteristic, – a compartment with a helium leak is simulated by two calculated volumes. In the domed space, the error in calculating helium concentration does not exceed 10%. The maximum error in determining the pressure (integral characteristic) before the sprinkler system starts functioning is ~ 5%.

The SIMCO code version 1.0 was developed at the Experimental Scientific-Research and Methodology Center “Simulation Systems” (SSL). Based on the results of testing and verification calculations, it may be concluded that physico-mathematical code models quite adequately describe the processes of heat/mass transfer in the containment. The described SIMCO code version can be used to analyze the totality of thermophysical and physico-chemical processes in nuclear power facilities with containments.