Development of unstructured mesh-based numerical method for sodium–water reaction phenomenon in steam generators of sodium-cooled fast reactors

ABSTRACT

When pressurized water or vapor leaks from a failed heat transfer tube in a steam generator of sodium-cooled fast reactors, a high-velocity and high-temperature jet with sodium–water chemical reaction may cause wastage on the adjacent tubes. For safety assessment of the steam generator, a computational fluid dynamics code called SERAPHIM calculating compressible multicomponent multiphase flow with sodium–water chemical reaction has been developed. The original SERAPHIM code is based on the finite difference method. In this study, unstructured mesh-based numerical method for the SERAPHIM code was developed to advance a numerical accuracy for the complex-shaped domain including multiple heat transfer tubes. Numerical analysis of an underexpanded jet experiment was performed as part of validation of the unstructured mesh-based SERAPHIM code. The calculated pressure profile showed good agreement with the experimental data. To investigate the effect of the introduction of the unstructured mesh and to confirm applicability of the numerical method for the actual situation, water vapor discharging into liquid sodium was analyzed. The calculated behavior of the reacting jet agreed with the previous experimental knowledge. It was demonstrated that the proposed numerical method could be applicable to evaluation of the sodium–water reaction phenomenon.

KEYWORDS:

• Sodium-cooled fast reactor
• steam generator
• multiphysics
• sodium–water reaction
• multiphase flow
• computational fluid dynamics
• numerical analysis
• unstructured mesh

1. Introduction

In a steam generator (SG) of sodium-cooled fast reactors, heat exchange takes place between liquid sodium in a shell side and water inside a heat transfer tube. When pressurized water or vapor leaks from a failed heat transfer tube, a high-velocity, high-temperature, and corrosive jet with sodium–water chemical reaction is formed in the shell side as shown in Figure 1. It is known that the reacting jet may cause wastage on the adjacent tubes. Wastage is due to erosion, flow accelerated corrosion, or combination of them. Significant progress of wastage or degradation of a mechanical strength by temperature rise may result in secondary failure (failure propagation) on the adjacent tubes. The tube failure propagation will lead to expansion of the damage and long-term shutdown of the plant.

Figure 1. Reacting jet under tube failure accident.

Prevention of the failure propagation is one of the major concerns in designing the SG. Our goal is to construct a computational method evaluating possibility of failure propagation occurrence. A computational method has a great advantage over the conventional evaluation method based on exhaustive demonstration experiments in terms of a cost, development period, and range of application. As shown in Figure 2, a multiphysics computational evaluation system comprising the multiple computational fluid dynamics (CFD) codes was newly constructed for design and safety evaluation of the SG in our study. The evaluation system includes modeling of multiphase flow, heat and mass transfer, chemical reaction, stress generation, variation of material properties, and failure of material. As one component of the system, a CFD code called SERAPHIM (Sodium-watEr Reaction Analysis: PHysics of Interdisciplinary Multi-phase flow) was developed to calculate compressible multicomponent multiphase flow with sodium–water chemical reaction [1–6]. The SERAPHIM code is based on the mechanistic numerical models. The multidimensional multifluid model considering compressibility and a numerical model for chemical reaction which occurs at the interface between water vapor and liquid sodium was developed in our previous studies. This enabled us to calculate the reacting multiphase flow with water, liquid sodium, and multicomponent gas.

Figure 2. Multiphysics evaluation system for tube failure accident.

The original SERAPHIM code uses the finite difference method. When a structured mesh is applied to the domain in which the heat transfer tube bundle exists, the tubes are modeled in a staircase pattern. In this case, significantly fine meshing is required to achieve enough resolution for the tubes, leading to large number of computational cells. On the other hand, an unstructured mesh can represent the curvilinear surface of the tubes with the relatively low number of the computational cells. In addition, the structured mesh has a problem to analyze the situation including multiple jets after the secondary failure. The secondary leakage will be formed diagonally to the axis of the primary jet. In this case, correct meshing for the secondary leakage on the surface of the target tube becomes extremely difficult. From these reasons, an unstructured mesh-based numerical method for the compressible multicomponent multiphase flow with sodium–water chemical reaction was developed and introduced into the SERAPHIM code in this study.

This paper focuses on construction of the unstructured mesh-based numerical method, investigation into fundamental validity of the method through the analysis of an underexpanded jet, and the test analysis on the reacting jet. Water vapor blowing out from an opening of the failed tube forms the underexpanded jet because a critical pressure of water vapor inside the tube is higher than a sodium-side ambient pressure. The gaseous underexpanded jet is one of the key phenomena under the tube failure accident. Numerical analysis of a three-dimensional underexpanded air jet in air was first performed to investigate validity of the unstructured mesh-based governing equations and the solution method considering compressibility. Kuehner et al. [7] and Woodmansee et al. [8] measured the profile of the pressure along the centerline of the underexpanded air jet in air by using the coherent anti-Stokes Raman scattering technique. The calculated pressure was compared to their experimental data. In addition, the test analysis on water vapor discharging into liquid sodium was performed to confirm a basic capability of the numerical method and to compare the effect of two types of meshing approaches. The hypothetical computational domain constructed for this test analysis consists of the simulated horizontal heat transfer tube (target tube) and the discharging nozzle below the target tube. The domain was initially filled with liquid sodium. Water vapor goes into the sodium pool vertically from the discharging nozzle and reacts with liquid sodium. The calculated behavior of the reacting jet is discussed in the following section.

2. Numerical methods

2.1. Governing equations

The multiphase flow under the tube failure accident is calculated by the multidimensional multifluid model considering compressibility. The first phase is a multicomponent gas including water vapor, sodium hydroxide (NaOH) aerosol, hydrogen, sodium vapor, and NaOH vapor. The NaOH aerosol and the hydrogen are the reaction products as described in Section 2.3 and are considered as one component of the gas phase. Evaporation of the NaOH aerosol by the reaction heat is calculated and the generated NaOH vapor is also considered as one component. The second phase is liquid sodium. The governing equations for the unstructured mesh are formulated as follows:

• Equation of mass conservation (1) $$\begin{array}{ccc}& & \frac{\partial }{\partial t}{\int }_{\Omega }{\alpha }_g{\rho }_{g}d\Omega +{\int }_S{\alpha }_g{\rho }_g{\mathbf{u}}_g·\mathbf{n}dS\hfill \\ & & ={\int }_{\Omega }\left({\Gamma }^e-{\Gamma }^c+G^{\mathrm{sf}}+G^{\mathrm{dif}}\right)d\Omega \hfill \end{array}$$(1) (2) $$\begin{array}{ccc}& & \frac{\partial }{\partial t}{\int }_{\Omega }{\alpha }_l{\rho }_{l}d\Omega +{\int }_S{\alpha }_l{\rho }_l{\mathbf{u}}_l·\mathbf{n}dS\hfill \\ & & ={\int }_{\Omega }\left(-{\Gamma }^e+{\Gamma }^c-G^{\mathrm{sf}}\right)d\Omega \hfill \end{array}$$(2)

• Equation of momentum conservation (3) $$\begin{array}{ccc}& & \frac{\partial }{\partial t}{\int }_{\Omega }{\alpha }_g{\rho }_g{\mathbf{u}}_{g}d\Omega +{\int }_S{\alpha }_g{\rho }_g{\mathbf{u}}_g{\mathbf{u}}_g·\mathbf{n}dS\hfill \\ & & =-{\int }_S{\alpha }_{g}p\mathbf{i}·\mathbf{n}dS+{\int }_{\Omega }\nabla ·\left({\alpha }_g{\tau }_g\right)d\Omega \hfill \\ & & +{\int }_{\Omega }({\alpha }_g{\rho }_g\mathbf{g}+{\mathbf{F}}^d+{\Gamma }^e{\mathbf{u}}_l-{\Gamma }^c{\mathbf{u}}_g+G^{\mathrm{sf}}{\mathbf{u}}_l\hfill \\ & & -G^{\mathrm{dif}}{\mathbf{u}}_g)d\Omega \hfill \end{array}$$(3) (4) $$\begin{array}{ccc}& & \frac{\partial }{\partial t}{\int }_{\Omega }{\alpha }_l{\rho }_l{\mathbf{u}}_{l}d\Omega +{\int }_S{\alpha }_l{\rho }_l{\mathbf{u}}_l{\mathbf{u}}_l·\mathbf{n}dS\hfill \\ & & =-{\int }_S{\alpha }_{l}p\mathbf{i}·\mathbf{n}dS+{\int }_{\Omega }\nabla ·\left({\alpha }_l{\tau }_l\right)d\Omega \hfill \\ & & +{\int }_{\Omega }\left({\alpha }_l{\rho }_l\mathbf{g}-{\mathbf{F}}^d-{\Gamma }^e{\mathbf{u}}_l+{\Gamma }^c{\mathbf{u}}_g-G^{\mathrm{sf}}{\mathbf{u}}_l\right)d\Omega \hfill \end{array}$$(4)

• Equation of energy conservation (5) $$\begin{array}{ccc}& & \frac{\partial }{\partial t}{\int }_{\Omega }{\alpha }_g{\rho }_g{h}_{g}d\Omega +{\int }_S{\alpha }_g{\rho }_g{h}_g{\mathbf{u}}_g·\mathbf{n}dS\hfill \\ & & ={\int }_S{\alpha }_g{\lambda }_g\nabla T_g·\mathbf{n}dS+{\int }_{\Omega }\{{\alpha }_g\frac{Dp}{Dt}+{\alpha }_g{\Phi }_g\hfill \\ & & +aH\left(T_l-T_g\right)+{\Gamma }^e{h}_{g\left(Na\right)}-{\Gamma }^c\left(h_g-i_l\right)\hfill \\ & & +Q^{\mathrm{sf}}+Q^{\mathrm{dif}}\}d\Omega \hfill \end{array}$$(5) (6) $$\begin{array}{ccc}& & \frac{\partial }{\partial t}{\int }_{\Omega }{\alpha }_l{\rho }_l{h}_{l}d\Omega +{\int }_S{\alpha }_l{\rho }_l{h}_l{\mathbf{u}}_l·\mathbf{n}dS\hfill \\ & & ={\int }_S{\alpha }_l{\lambda }_l\nabla T_l·\mathbf{n}dS+{\int }_{\Omega }\{{\alpha }_l\frac{Dp}{Dt}+{\alpha }_l{\Phi }_l\hfill \\ & & +aH\left(T_g-T_l\right)-{\Gamma }^e\left(h_l+i_l\right)+{\Gamma }^c{h}_l-G^{\mathrm{sf}}{h}_l\}d\Omega \hfill \end{array}$$(6)

where α is the volume fraction, ρ the density, u the velocity vector, Γ^e the evaporation rate, Γ^c the condensation rate, G^sf the mass generation rate by surface reaction, G^dif the source term by species diffusion, p the pressure, τ the viscous stress tensor, g the gravity vector, F^d the interfacial drag force, h the enthalpy, Φ the dissipation function, λ the thermal conductivity, T the temperature, a the interfacial area density, H the coefficient of heat transfer, i the latent heat, h_g(Na) the enthalpy of vapor sodium, Q^sf the heat generation rate by surface reaction, Q^dif the source term by species diffusion, Ω the volume, and S the surface area. Subscript g stands for the gas phase and l the liquid phase. The source terms G^dif and Q^dif are written as (7) $$\begin{array}{c}\hfill G^{\mathrm{dif}}=\sum _j\left(\nabla ·{\alpha }_g{\rho }_g{D}_{mj}\nabla Y_j\right)\end{array}$$(7) (8) $$\begin{array}{c}\hfill Q^{\mathrm{dif}}=\sum _j\left(\nabla ·{\alpha }_g{\rho }_g{D}_{mj}{h}_j\nabla Y_j\right)\end{array}$$(8) where j is the component of the gas, D_m the effective coefficient of diffusion, and Y the mass fraction. Equation (7)(7) $$\begin{array}{c}\hfill G^{\mathrm{dif}}=\sum _j\left(\nabla ·{\alpha }_g{\rho }_g{D}_{mj}\nabla Y_j\right)\end{array}$$(7) corresponds to the first term on the right-hand side in Equation (9)(9) $$\begin{array}{ccc}& & \frac{\partial }{\partial t}{\int }_{\Omega }{\alpha }_g{\rho }_g{Y}_{j}d\Omega +{\int }_S{\alpha }_g{\rho }_g{Y}_j{\mathbf{u}}_g·\mathbf{n}dS\hfill \\ & & ={\int }_S{\alpha }_g{\rho }_g{D}_{mj}\nabla Y_j·\mathbf{n}dS+{\int }_{\Omega }\left({\Gamma }_j^e-{\Gamma }_j^c+{\gamma }_j^{\mathrm{sf}}\right)d\Omega \hfill \end{array}$$(9) to be described. Transport of the species j in the multicomponent gas is calculated from the advection–diffusion equation: (9) $$\begin{array}{ccc}& & \frac{\partial }{\partial t}{\int }_{\Omega }{\alpha }_g{\rho }_g{Y}_{j}d\Omega +{\int }_S{\alpha }_g{\rho }_g{Y}_j{\mathbf{u}}_g·\mathbf{n}dS\hfill \\ & & ={\int }_S{\alpha }_g{\rho }_g{D}_{mj}\nabla Y_j·\mathbf{n}dS+{\int }_{\Omega }\left({\Gamma }_j^e-{\Gamma }_j^c+{\gamma }_j^{\mathrm{sf}}\right)d\Omega \hfill \end{array}$$(9) where γ^sf is the source term by surface reaction.

The equation of state is considered in addition to Equations (1)(1) $$\begin{array}{ccc}& & \frac{\partial }{\partial t}{\int }_{\Omega }{\alpha }_g{\rho }_{g}d\Omega +{\int }_S{\alpha }_g{\rho }_g{\mathbf{u}}_g·\mathbf{n}dS\hfill \\ & & ={\int }_{\Omega }\left({\Gamma }^e-{\Gamma }^c+G^{\mathrm{sf}}+G^{\mathrm{dif}}\right)d\Omega \hfill \end{array}$$(1) –(9(9) $$\begin{array}{ccc}& & \frac{\partial }{\partial t}{\int }_{\Omega }{\alpha }_g{\rho }_g{Y}_{j}d\Omega +{\int }_S{\alpha }_g{\rho }_g{Y}_j{\mathbf{u}}_g·\mathbf{n}dS\hfill \\ & & ={\int }_S{\alpha }_g{\rho }_g{D}_{mj}\nabla Y_j·\mathbf{n}dS+{\int }_{\Omega }\left({\Gamma }_j^e-{\Gamma }_j^c+{\gamma }_j^{\mathrm{sf}}\right)d\Omega \hfill \end{array}$$(9) ). (10) $$\begin{array}{c}\hfill p=z\frac{{\rho }_g}{M_W}R{T}_g\end{array}$$(10) where z is the coefficient of compressibility, M_W the molecular weight, and R the universal gas constant. The ideal gas approximation (z = 1) gives certain errors for the high-pressure gas such as the water vapor in the SG. In our method, the coefficient of compressibility is calculated from the modified Benedict–Webb–Rubin equation [9] to obtain the actual gas-phase density.

The interfacial drag force term and the interfacial area density in the governing equations are estimated by the model proposed by Autruffe et al. [3,10] and Nigmatulin et al. [3,11], respectively. The Nigmatulin model includes empirical parameter representing the number density of the bubbles or the liquid droplet. In our previous researches, the parameter was decided so as to reproduce the temperature field formed by the reacting jet from sensitivity analysis. In other words, the waviness effect appearing in the actual reacting jet is considered in this empirical parameter.

In the case of applying the SERAPHIM code to the actual situation, the above governing equations are solved for the region in which the heat transfer tube bundle exists. In this case, the reacting jet impinges on the multiple tubes and spreads vigorously. No turbulence model is taken into account in the governing equations because the effect of spreading by impingement on the tubes seems to be larger than the turbulence diffusion.

The time-derivative terms in Equations (1)(1) $$\begin{array}{ccc}& & \frac{\partial }{\partial t}{\int }_{\Omega }{\alpha }_g{\rho }_{g}d\Omega +{\int }_S{\alpha }_g{\rho }_g{\mathbf{u}}_g·\mathbf{n}dS\hfill \\ & & ={\int }_{\Omega }\left({\Gamma }^e-{\Gamma }^c+G^{\mathrm{sf}}+G^{\mathrm{dif}}\right)d\Omega \hfill \end{array}$$(1) –(6(6) $$\begin{array}{ccc}& & \frac{\partial }{\partial t}{\int }_{\Omega }{\alpha }_l{\rho }_l{h}_{l}d\Omega +{\int }_S{\alpha }_l{\rho }_l{h}_l{\mathbf{u}}_l·\mathbf{n}dS\hfill \\ & & ={\int }_S{\alpha }_l{\lambda }_l\nabla T_l·\mathbf{n}dS+{\int }_{\Omega }\{{\alpha }_l\frac{Dp}{Dt}+{\alpha }_l{\Phi }_l\hfill \\ & & +aH\left(T_g-T_l\right)-{\Gamma }^e\left(h_l+i_l\right)+{\Gamma }^c{h}_l-G^{\mathrm{sf}}{h}_l\}d\Omega \hfill \end{array}$$(6) ) and (9(9) $$\begin{array}{ccc}& & \frac{\partial }{\partial t}{\int }_{\Omega }{\alpha }_g{\rho }_g{Y}_{j}d\Omega +{\int }_S{\alpha }_g{\rho }_g{Y}_j{\mathbf{u}}_g·\mathbf{n}dS\hfill \\ & & ={\int }_S{\alpha }_g{\rho }_g{D}_{mj}\nabla Y_j·\mathbf{n}dS+{\int }_{\Omega }\left({\Gamma }_j^e-{\Gamma }_j^c+{\gamma }_j^{\mathrm{sf}}\right)d\Omega \hfill \end{array}$$(9) ) are approximated by the forward difference scheme. The second-order total variation diminishing (TVD) scheme is an effective measure to reproduce the behavior of an underexpanded jet. Darwish and Moukalled [12] proposed the second-order TVD scheme for an unstructured mesh. We applied their method to the convective terms in the equation of momentum conservation. In our present work, application of the second-order TVD scheme for the convective terms in other equations caused numerical instability. Therefore, the convective terms in other equations were discretized by using the first-order upwind difference scheme at this stage. The second-order central difference scheme was applied to the viscous, heat conduction and diffusion terms.

2.2. HSMAC method for compressible multiphase flow

The governing equations are solved by the HSMAC (Highly Simplified Marker And Cell) method modified for compressible multiphase flows [2,4]. The HSMAC method for incompressible single-phase flow was developed by Hirt and Cook [13]. In this method, the pressure and the velocity are corrected so as to satisfy continuity by using the Newton–Raphson method. The HSMAC method is very useful for parallel computation because there is no need to solve a simultaneous linear equation. Formulation of the HSMAC method for compressible multiphase flows under an unstructured mesh is described below. The recurrence equation for the pressure correction is written as (11) $$\begin{array}{c}\hfill p^{m+1}=p^m+\delta p=p^m-{D^m/\phantom{D^m\left(\frac{\partial D}{\partial p}\right)}\left(\frac{\partial D}{\partial p}\right)}^m\end{array}$$(11) where superscript m stands for the number of iterations. The variable D approaches to 0 from iterative calculation. In the case of multiphase flows, similarly to the method proposed by Matsumoto and Murai [14], we define the variable D by (12) $$\begin{array}{c}\hfill D=-\sum _k\frac{\partial {\alpha }_k}{\partial t}\end{array}$$(12)

The equations of mass conservation (1(1) $$\begin{array}{ccc}& & \frac{\partial }{\partial t}{\int }_{\Omega }{\alpha }_g{\rho }_{g}d\Omega +{\int }_S{\alpha }_g{\rho }_g{\mathbf{u}}_g·\mathbf{n}dS\hfill \\ & & ={\int }_{\Omega }\left({\Gamma }^e-{\Gamma }^c+G^{\mathrm{sf}}+G^{\mathrm{dif}}\right)d\Omega \hfill \end{array}$$(1) ) and (2(2) $$\begin{array}{ccc}& & \frac{\partial }{\partial t}{\int }_{\Omega }{\alpha }_l{\rho }_{l}d\Omega +{\int }_S{\alpha }_l{\rho }_l{\mathbf{u}}_l·\mathbf{n}dS\hfill \\ & & ={\int }_{\Omega }\left(-{\Gamma }^e+{\Gamma }^c-G^{\mathrm{sf}}\right)d\Omega \hfill \end{array}$$(2) ) are rewritten as (13) $$\begin{array}{c}\hfill \frac{\partial \left({\alpha }_k{\rho }_k\right)}{\partial t}\Delta \Omega +\sum _c\left\{⟨{\alpha }_k{\rho }_k⟩{\mathbf{u}}_k·\mathbf{n}{S}_c\right\}={\Gamma }_k^{\mathrm{sum}}\Delta \Omega \end{array}$$(13) (14) $$\begin{array}{c}\hfill {\rho }_k\frac{\partial {\alpha }_k}{\partial t}+{\alpha }_k\frac{\partial {\rho }_k}{\partial t}+\frac{1}{\Delta \Omega }\sum _c\left\{⟨{\alpha }_k{\rho }_k⟩{\mathbf{u}}_k·\mathbf{n}{S}_c\right\}={\Gamma }_k^{\mathrm{sum}}\end{array}$$(14) where ΔΩ is the volume of a cell and S_c is the area of the one surface of the cell. The angle bracket stands for the mean value in each surface of the cell. Substituting Equation (14)(14) $$\begin{array}{c}\hfill {\rho }_k\frac{\partial {\alpha }_k}{\partial t}+{\alpha }_k\frac{\partial {\rho }_k}{\partial t}+\frac{1}{\Delta \Omega }\sum _c\left\{⟨{\alpha }_k{\rho }_k⟩{\mathbf{u}}_k·\mathbf{n}{S}_c\right\}={\Gamma }_k^{\mathrm{sum}}\end{array}$$(14) into Equation (12)(12) $$\begin{array}{c}\hfill D=-\sum _k\frac{\partial {\alpha }_k}{\partial t}\end{array}$$(12) , D is rewritten as (15) $$\begin{array}{ccc}\hfill D& =& \sum _k[\frac{{\alpha }_k}{{\rho }_k}\frac{\partial {\rho }_k}{\partial t}+\frac{1}{{\rho }_k\Delta \Omega }\sum _c\left\{⟨{\alpha }_k{\rho }_k⟩{\mathbf{u}}_k·\mathbf{n}{S}_c\right\}\hfill \\ & & -\frac{{\Gamma }_k^{\mathrm{sum}}}{{\rho }_k}]\hfill \end{array}$$(15) ∂ρ_k/∂t in Equation (13)(13) $$\begin{array}{c}\hfill \frac{\partial \left({\alpha }_k{\rho }_k\right)}{\partial t}\Delta \Omega +\sum _c\left\{⟨{\alpha }_k{\rho }_k⟩{\mathbf{u}}_k·\mathbf{n}{S}_c\right\}={\Gamma }_k^{\mathrm{sum}}\Delta \Omega \end{array}$$(13) can be estimated from Equation (8)(8) $$\begin{array}{c}\hfill Q^{\mathrm{dif}}=\sum _j\left(\nabla ·{\alpha }_g{\rho }_g{D}_{mj}{h}_j\nabla Y_j\right)\end{array}$$(8) as follows: (16) $$\begin{array}{c}\hfill \frac{\partial {\rho }_k}{\partial t}=\frac{\partial {\rho }_k}{\partial p}\frac{\partial p}{\partial t}+\frac{\partial {\rho }_k}{\partial T_k}\frac{\partial T_k}{\partial t}+\frac{\partial {\rho }_k}{\partial M_W}\frac{\partial M_W}{\partial t}\end{array}$$(16) Variation of D in the iterative calculation is then written as (17) $$\begin{array}{ccc}\hfill \delta D& =& D^{n+1}-D^{*}=\sum _k\frac{{\alpha }_k^n}{{\rho }_k^n}\frac{\partial {\rho }_k}{\partial p}\frac{\delta p}{\Delta t}+\sum _k\frac{{\alpha }_k^n}{{\rho }_k^n}\frac{\partial {\rho }_k}{\partial T_k}\frac{\delta T_k}{\Delta t}\hfill \\ & & +\frac{{\alpha }_g^n}{{\rho }_g^n}\frac{\partial {\rho }_g}{\partial M_W}\frac{\delta M_W}{\Delta t}+\sum _k[\frac{1}{{\rho }_k^n\Delta \Omega }\hfill \\ & & \times \sum _c\left\{⟨{\alpha }_k^n{\rho }_k^n⟩\delta {\mathbf{u}}_k^{n+1}·\mathbf{n}{S}_c\right\}]\hfill \end{array}$$(17) where (18) $$\begin{array}{c}\hfill \delta p=p^{n+1}-p^{*}\end{array}$$(18) (19) $$\begin{array}{c}\hfill \delta T_k=T_k^{n+1}-T_k^{*}\end{array}$$(19) (20) $$\begin{array}{c}\hfill \delta M_W=M_W^{n+1}-M_W^{*}\end{array}$$(20) (21) $$\begin{array}{c}\hfill \delta {\mathbf{u}}_k={\mathbf{u}}_k^{n+1}-{\mathbf{u}}_k^{*}\end{array}$$(21) and superscript n + 1 stands for the current time step and ^* the predictive value in the iteration. When the solution converges, D becomes 0 (Dⁿ⁺¹ = 0). This means that Equation (17)(17) $$\begin{array}{ccc}\hfill \delta D& =& D^{n+1}-D^{*}=\sum _k\frac{{\alpha }_k^n}{{\rho }_k^n}\frac{\partial {\rho }_k}{\partial p}\frac{\delta p}{\Delta t}+\sum _k\frac{{\alpha }_k^n}{{\rho }_k^n}\frac{\partial {\rho }_k}{\partial T_k}\frac{\delta T_k}{\Delta t}\hfill \\ & & +\frac{{\alpha }_g^n}{{\rho }_g^n}\frac{\partial {\rho }_g}{\partial M_W}\frac{\delta M_W}{\Delta t}+\sum _k[\frac{1}{{\rho }_k^n\Delta \Omega }\hfill \\ & & \times \sum _c\left\{⟨{\alpha }_k^n{\rho }_k^n⟩\delta {\mathbf{u}}_k^{n+1}·\mathbf{n}{S}_c\right\}]\hfill \end{array}$$(17) represents the relationship between the divergence D^* and the pressure correction δp. From Equation (17)(17) $$\begin{array}{ccc}\hfill \delta D& =& D^{n+1}-D^{*}=\sum _k\frac{{\alpha }_k^n}{{\rho }_k^n}\frac{\partial {\rho }_k}{\partial p}\frac{\delta p}{\Delta t}+\sum _k\frac{{\alpha }_k^n}{{\rho }_k^n}\frac{\partial {\rho }_k}{\partial T_k}\frac{\delta T_k}{\Delta t}\hfill \\ & & +\frac{{\alpha }_g^n}{{\rho }_g^n}\frac{\partial {\rho }_g}{\partial M_W}\frac{\delta M_W}{\Delta t}+\sum _k[\frac{1}{{\rho }_k^n\Delta \Omega }\hfill \\ & & \times \sum _c\left\{⟨{\alpha }_k^n{\rho }_k^n⟩\delta {\mathbf{u}}_k^{n+1}·\mathbf{n}{S}_c\right\}]\hfill \end{array}$$(17) , we obtain a final form of the pressure correction: (22) $$\begin{array}{c}\hfill \delta p=-\omega \frac{D^{*}+{\sum }_k\left\{\frac{{\alpha }_k^n}{{\rho }_k^n}\left(\frac{\partial {\rho }_k}{\partial T_k}\frac{\delta T_k}{\Delta t}+\frac{\partial {\rho }_k}{\partial M}\frac{\delta M}{\Delta t}\right)\right\}}{{\sum }_k\left[\frac{{\alpha }_k^n}{{\rho }_k^n}\left\{\frac{\partial {\rho }_k}{\partial p}\frac{1}{\Delta t}+\frac{1}{{\alpha }_k^n\Delta \Omega }{\sum }_c\left\{⟨{\alpha }_k^n{\rho }_k^n⟩\delta {\mathbf{u}}_k·\mathbf{n}{S}_c\right\}\right\}\right]}\end{array}$$(22) where ω is a relaxation factor. The T^{n + 1}_k and M^{n + 1}_W in Equations (19)(19) $$\begin{array}{c}\hfill \delta T_k=T_k^{n+1}-T_k^{*}\end{array}$$(19) and (20(20) $$\begin{array}{c}\hfill \delta M_W=M_W^{n+1}-M_W^{*}\end{array}$$(20) ) are estimated by latest values in the iteration. From the equations of momentum conservation (3(3) $$\begin{array}{ccc}& & \frac{\partial }{\partial t}{\int }_{\Omega }{\alpha }_g{\rho }_g{\mathbf{u}}_{g}d\Omega +{\int }_S{\alpha }_g{\rho }_g{\mathbf{u}}_g{\mathbf{u}}_g·\mathbf{n}dS\hfill \\ & & =-{\int }_S{\alpha }_{g}p\mathbf{i}·\mathbf{n}dS+{\int }_{\Omega }\nabla ·\left({\alpha }_g{\tau }_g\right)d\Omega \hfill \\ & & +{\int }_{\Omega }({\alpha }_g{\rho }_g\mathbf{g}+{\mathbf{F}}^d+{\Gamma }^e{\mathbf{u}}_l-{\Gamma }^c{\mathbf{u}}_g+G^{\mathrm{sf}}{\mathbf{u}}_l\hfill \\ & & -G^{\mathrm{dif}}{\mathbf{u}}_g)d\Omega \hfill \end{array}$$(3) ) and (4(4) $$\begin{array}{ccc}& & \frac{\partial }{\partial t}{\int }_{\Omega }{\alpha }_l{\rho }_l{\mathbf{u}}_{l}d\Omega +{\int }_S{\alpha }_l{\rho }_l{\mathbf{u}}_l{\mathbf{u}}_l·\mathbf{n}dS\hfill \\ & & =-{\int }_S{\alpha }_{l}p\mathbf{i}·\mathbf{n}dS+{\int }_{\Omega }\nabla ·\left({\alpha }_l{\tau }_l\right)d\Omega \hfill \\ & & +{\int }_{\Omega }\left({\alpha }_l{\rho }_l\mathbf{g}-{\mathbf{F}}^d-{\Gamma }^e{\mathbf{u}}_l+{\Gamma }^c{\mathbf{u}}_g-G^{\mathrm{sf}}{\mathbf{u}}_l\right)d\Omega \hfill \end{array}$$(4) ), δu_k in Equation (21)(21) $$\begin{array}{c}\hfill \delta {\mathbf{u}}_k={\mathbf{u}}_k^{n+1}-{\mathbf{u}}_k^{*}\end{array}$$(21) is given by (23) $$\begin{array}{ccc}\hfill \delta {\mathbf{u}}_k& =& {\mathbf{u}}_k^{n+1}-{\mathbf{u}}_k^{*}\hfill \\ & =& \frac{\Delta t}{{\rho }_k\Delta \Omega }\sum _c\left\{⟨{\alpha }_k^n\delta p⟩\mathbf{i}·\mathbf{n}{S}_c\right\}\hfill \end{array}$$(23)

The m + 1th pressure is calculated from Equations (11)(11) $$\begin{array}{c}\hfill p^{m+1}=p^m+\delta p=p^m-{D^m/\phantom{D^m\left(\frac{\partial D}{\partial p}\right)}\left(\frac{\partial D}{\partial p}\right)}^m\end{array}$$(11) , (22(22) $$\begin{array}{c}\hfill \delta p=-\omega \frac{D^{*}+{\sum }_k\left\{\frac{{\alpha }_k^n}{{\rho }_k^n}\left(\frac{\partial {\rho }_k}{\partial T_k}\frac{\delta T_k}{\Delta t}+\frac{\partial {\rho }_k}{\partial M}\frac{\delta M}{\Delta t}\right)\right\}}{{\sum }_k\left[\frac{{\alpha }_k^n}{{\rho }_k^n}\left\{\frac{\partial {\rho }_k}{\partial p}\frac{1}{\Delta t}+\frac{1}{{\alpha }_k^n\Delta \Omega }{\sum }_c\left\{⟨{\alpha }_k^n{\rho }_k^n⟩\delta {\mathbf{u}}_k·\mathbf{n}{S}_c\right\}\right\}\right]}\end{array}$$(22) ), and (23(23) $$\begin{array}{ccc}\hfill \delta {\mathbf{u}}_k& =& {\mathbf{u}}_k^{n+1}-{\mathbf{u}}_k^{*}\hfill \\ & =& \frac{\Delta t}{{\rho }_k\Delta \Omega }\sum _c\left\{⟨{\alpha }_k^n\delta p⟩\mathbf{i}·\mathbf{n}{S}_c\right\}\hfill \end{array}$$(23) ). The velocity is updated after the pressure correction. When the residual calculated from Equation (12)(12) $$\begin{array}{c}\hfill D=-\sum _k\frac{\partial {\alpha }_k}{\partial t}\end{array}$$(12) converges, we move to the calculation in the next time step.

2.3. Chemical reaction model

The mass generation rate by the chemical reaction at the interface between water vapor and liquid sodium is calculated by a surface reaction model developed in our previous studies. The surface reaction model is outlined in this section. The following reaction is dominant between water vapor and liquid sodium [1,2]. (24) $$\begin{array}{c}\hfill \mathrm{Na}+{\mathrm{H}}_2\mathrm{O}\to \mathrm{NaOH}+\frac{1}{2}{\mathrm{H}}_2\end{array}$$(24)

Figure 3 illustrates the surface reaction expressed by Equation (24)(24) $$\begin{array}{c}\hfill \mathrm{Na}+{\mathrm{H}}_2\mathrm{O}\to \mathrm{NaOH}+\frac{1}{2}{\mathrm{H}}_2\end{array}$$(24) . The surface reaction model is based on the assumption of infinite reaction rate. In other words, progress of the chemical reaction at the gas–liquid interface is limited by the mass flux of the reactant gas toward the interface. The mass flux of the reactant gas j is written as (25) $$\begin{array}{c}\hfill {\gamma }_j^{\mathrm{sf}}=\mathrm{Sh}\frac{D_{\mathrm{mj}}}{l}{\rho }_g{Y}_{j}a\end{array}$$(25) where Sh is the Sherwood number and l the characteristic length. An analogy between heat and mass transfer is applied to Equation (25)(25) $$\begin{array}{c}\hfill {\gamma }_j^{\mathrm{sf}}=\mathrm{Sh}\frac{D_{\mathrm{mj}}}{l}{\rho }_g{Y}_{j}a\end{array}$$(25) to replace some unknown parameters included in this equation. Equation (25)(25) $$\begin{array}{c}\hfill {\gamma }_j^{\mathrm{sf}}=\mathrm{Sh}\frac{D_{\mathrm{mj}}}{l}{\rho }_g{Y}_{j}a\end{array}$$(25) is rewritten as (26) $$\begin{array}{c}\hfill {\gamma }_j^{\mathrm{sf}}=-\mathrm{L}{\mathrm{e}}^{b-1}\frac{H_{\mathrm{gl}}}{C_{\mathrm{pg}}}{Y}_{j}a\end{array}$$(26) where Le is the Lewis number, b the empirical constant, and C_pg the gas-phase constant pressure specific heat. The method to calculate each parameter can be found in our previous reports [6].

Figure 3. Chemical reaction at interface between water vapor and liquid sodium.

3. Numerical analysis

3.1. Analysis of underexpanded air jet in air

Water vapor blowing out from an opening of the failed tube forms an underexpanded jet because a critical pressure of water vapor inside the tube is higher than a sodium-side ambient pressure. Figure 4 shows the well-known structure of the underexpanded jet. The underexpanded jet is characterized by the appearance of a shock wave. A Mach disk, which is one of the shock wave, is almost normal to a centerline of the jet. It is known that the various physical quantities vary sharply at the Mach disk. The experiment on the underexpanded air jet in air by Kuehner et al. [7] and Woodmansee et al. [8] was analyzed to confirm validity of the above-mentioned governing equations and solution method considering compressibility. We considered the three-dimensional cylindrical domain shown in Figure 5(a). The domain is filled with air at the pressure of 0.1 MPa and the temperature of 23.85 °C. Pressurized air goes into the region from the inflow boundary corresponding to nozzle exit. The stagnation pressure p₀ and the temperature T₀ of pressurized air was set to 0.617 MPa and 23.85 °C. Air at the inflow boundary is in the critical state. Assuming an isentropic flow, the physical quantities under the critical state are given by the following relationships: (27) $$u^{*}=\sqrt{\frac{2\gamma }{\gamma +1}\frac{p_0}{{\rho }_0}}$$(27) (28) $$p^{*}=p_0{\left(\frac{2}{\gamma +1}\right)}^{\frac{\gamma }{\gamma -1}}$$(28) (29) $$T^{*}=T_0\frac{2}{\gamma +1}$$(29) (30) $${\rho }^{*}={\rho }_0{\left(\frac{2}{\gamma +1}\right)}^{\frac{1}{\gamma -1}}$$(30) where superscript ^* stands for the critical state and subscript 0 stands for the stagnation state. The velocity u^* becomes a sound speed at the inflow boundary. In the case of p₀ = 0.617 MPa, the critical pressure p^* is 0.325 MPa from Equation (28)(28) $$p^{*}=p_0{\left(\frac{2}{\gamma +1}\right)}^{\frac{\gamma }{\gamma -1}}$$(28) . The pressure and the temperature at the side surface were assumed to be constant at 0.1 MPa (ambient pressure) and 23.85 °C. An underexpanded jet appears in this case because the critical pressure is higher than the ambient pressure. The gradient of the pressure and the temperature at the outflow boundary was set to be 0. The symmetric condition was applied to the surface of x = 0 and y = 0. Figure 5(b) shows the analysis mesh at the surface of z = 0. In a direction of z, the domain was divided into the equally spaced mesh with a cell 0.15625 mm on a side.

Figure 4. Structure of underexpanded jet.

Figure 5. Computational domain and mesh for underexpanded air jet in air.

Figure 6 shows the time-averaged distributions of the velocity and the Mach number. It can be seen that the air expands after exiting the nozzle and the Mach disk appears in the jet. The velocity reaches up to a supersonic upstream of the Mach disk. There is a region of the low velocity downstream of the Mach disk. The numerical result reproduces the structure of the underexpanded jet shown in Figure 4. Figure 7 shows the profile of the time-averaged pressure along the centerline of the jet. The pressure decreases with underexpansion and increases sharply at the Mach disk. The calculated pressure shows good agreement with the experimental data by Kuehner et al. or by Woodmansee et al. It was demonstrated that the proposed numerical method could represent the behavior of the underexpanded jets, which is one of the most basic phenomena under tube failure accident.

Figure 6. Numerical results of underexpanded air jet in air.

Figure 7. Pressure profiles along centerline of underexpanded air jet.

3.2. Analysis of water vapor discharging into sodium pool

The test analysis on water vapor discharging into liquid sodium was performed to confirm a capability of the numerical method and to compare the effect of two types of meshing approaches. The effect of use of the unstructured mesh must be investigated quantitatively through a comparison with experimental results. The work on numerical analysis of sodium–water reaction experiment conducted in the domain including heat transfer tube bundle will be performed in our future researches. The hypothetical domain shown in Figure 8(a) was constructed for this test analysis. The domain consists of the simulated horizontal heat transfer tube (target tube) and the discharging nozzle below the target tube. The domain was initially filled with the liquid sodium at the temperature of 470 °C. There is an 8-mm square-shaped inflow boundary at the upper end of the discharging nozzle. The water vapor at 470 °C and 10 MPa goes into the pool vertically from the inflow boundary with 100 m/s and the reacting flow impinges on the target tube. The liquid sodium flows in a vertical direction at 5 m/s from the bottom surface of the domain. The gas and liquid phases flow out from the top surface. The pressure at the top surface was set to be constant at 0.2 MPa. The temperature and the pressure conditions shown here is approximately same as the actual conditions of the SG. Figure 8(b) shows the two different computational meshes applied to this analysis. The computational mesh A is constructed with rectangular cells. The surface of the target tube is modeled in a staircase pattern. On the other hand, the curvilinear surface of the target tube is simulated well in the computational mesh B constructed with triangular prism cells. The division number for each side of the rectangular domain conforms completely between the meshes A and B. The surface of the target tube is divided by roughly the same number of computational cells. In this case, the number of computational cells is 52,408 in the mesh A and 121,168 in the mesh B.

Figure 8. Computational domain for water vapor discharging into sodium pool.

Figure 9 shows variation of the distributions of the void fraction and the gas-phase velocity on a vertical plane. This figure compares the numerical results under the two meshes. The numerical results show that the gas-phase goes upward around the target tube. The target tube in a staircase pattern increases the effect as an obstacle in the flow compared to that in the unstructured mesh. Therefore, the gas-phase velocity near the tube surface is relatively low at 0.15 seconds in the mesh A. The void fraction at the upper side of the target tube is also different between the meshes A and B at 0.15 seconds. Figure 10 shows the instantaneous distributions of the gas-phase temperature, the volume fraction of water vapor, the volume fraction of hydrogen, and the mass concentration of NaOH at 0.15 seconds. The volume fraction of water vapor is almost 1 above the discharging nozzle in both the meshes A and B. In the same region, the gas-phase temperature is relatively low due to the existence of the unreacted water vapor. The water vapor is used in chemical reaction with its rising. Instead, the reaction products, hydrogen and NaOH, appear and become the main species of the gas phase. The high-temperature region is formed around the unreacted water vapor and above the target tube due to the chemical reaction. At the upper side of the target tube, the significant difference of the distributions of the chemical species appears between the meshes A and B. The flow of water vapor is almost completely separated by the target tube in the mesh A while the water vapor flows behind the target tube in the mesh B. From this, the gas-phase temperature at the upper side of the target tube in the mesh B is higher than that in the mesh A. The temperature distribution around the target tube is extremely important in predicting occurrence of the tube failure propagation since the material strength of the tube depends on the ambient fluid temperature. The use of the unstructured mesh is highly effective on the prediction. The gas-phase temperature increased up to about 1270 °C by the chemical reaction in the mesh B. Kurihara et al. [15] measured the temperature distribution in the experiment on discharging of the water vapor in the liquid sodium pool in which the heat transfer tube bundle exists and reported that the peak temperature increased over 1200 °C. The numerical results demonstrate that the developed SERAPHIM code has a capability to predict the peak temperature formed by the sodium–water chemical reaction. Numerical analysis of the experiment on sodium–water reaction phenomenon with multi-tube configuration by using the unstructured mesh-based SERAPHIM code will be performed to enhance its reliability in our future work.

Figure 9. Numerical results of water vapor discharging into sodium pool (variation of void fraction and gas-phase velocity).

Figure 10. Numerical results of water vapor discharging into sodium pool (gas-phase temperature and chemical species at 0.015 seconds).

4. Conclusion

The unstructured mesh-based numerical method was developed and introduced into the SERAPHIM code to advance a numerical accuracy for the complex-shaped domain including multiple heat transfer tubes. The validity of the numerical method was investigated through some numerical analyses. First, it was demonstrated that the proposed governing equations and the solution method considering compressibility could represent underexpanded jets, which is one of the most basic phenomena under tube failure accident. The numerical results reproduced a well-known structure of the underexpanded jet. The calculated pressure showed good agreement with the experimental data. In the numerical analysis of water vapor discharging into liquid sodium, the SERAPHIM code could calculate the temperature field formed by the multiphase flow with chemical reaction. The use of the unstructured mesh is highly effective in predicting the distributions of the fluid temperature around the target tube. The calculated peak temperature agreed with the previous experimental knowledge. The proposed numerical method has a capability to evaluate the sodium–water reaction phenomenon.

Disclosure statement

No potential conflict of interest was reported by the authors.