The development and validation of a natural circulation analysis code for marine reactors

ABSTRACT

A code PNCMC (Program for Natural Circulation under Motion Conditions) has been developed for natural circulation simulation of marine reactors. The code is based on one-dimensional two-fluid model in noninertial frame of reference. The body force term in the momentum equation is considered as a time dependent function, which consists of gravity and inertial force induced by three-dimensional ship motion. Staggered mesh, finite volume method, semi-implicit first order upwind scheme and Successive Over Relaxation (SOR) method are used to discretize and solve two-phase mass, momentum and energy equations. Single-phase natural circulation experiments under rolling condition performed in Institute of nuclear and new energy technology of Tsinghua University and two-phase natural circulation experiments under rolling condition performed by Tan and colleagues are used to validate PNCMC. The validation results indicate that PNCMC is capable to investigate the single-phase and two-phase natural circulation under rolling motion.

KEYWORDS:

• Natural circulation
• marine reactor
• motion condition
• thermal hydraulics

1. Introduction

Innovative small (less than 300 MWe) and medium (between 300 and 700 MWe) sized reactors are highly valued in both developing and developed countries. The Intensified efforts on small and medium reactors (SMRs) designs and research generally aim to provide increased benefits in the area of safety and security, as well as power generation flexibility. Apart from being designed as land-based nuclear power plants, small and medium sized light water reactors (LWRs) can be mounted on a barge as FNPPs (floating nuclear power plants) or as propulsion power. For example, the SMART 330 MW integral pressurized water reactor was designed for seawater desalination, small-scale power generation and ship propulsion. KLT-40s is an upgraded version of ice-breaker ship reactor KLT-40, which is designed for FNPP. So, operating in the ocean becomes an important development trend of small and medium sized LWRs. Most notably, natural circulation is widely used in the design of small and medium sized LWRs, such as ABV-6M, VK-300, NuScale and IMR [1]. Almost all the small and medium sized LWRs introduce passive residual heat removal systems, and the essence of these systems is also natural circulation.

However, during operation in the ocean, the sea wave together with wind induce the barge rolling, pitching and heaving, and these motions affect the reactor thermal hydraulics, especially natural circulation behavior [2]. The effects of ship motion focus on two principles. First, inclination decreases the height difference between the source (core) and sink (steam generator [SG]), and thus decreases the thermal driving head of natural circulation; Second, inertial forces imparted by the ship motion in addition to the gravity act upon the primary coolant, which together with flow resistance change the natural circulation flow rate. The direction and magnitude of inertial forces imparted by rolling depend on the position of rolling axis, rolling acceleration and reactor system arrangement in three-dimensional space. So the natural circulation behaviors under ocean conditions are much more complicated.

Hence considering the nuclear safety, a ton of research activities aiming at investigating natural circulation under motion condition have been performed over the years. Japan Atomic Energy Research Institute (JAERI) conducted experiments using the first Japanese experimental ship Mutsu in various sea conditions to investigate the behavior of the whole marine reactor system from May 1991 to February 1992 [3]. They measured the parameters of wave height, wind force, ship motion (acceleration) and reactor thermal hydraulics; the results indicated that the ship motions hardly affect the operation parameters of primary loop (forced circulation loop), but prominently affect the steam generator water level, evaporator steam void fraction and recirculation flow rate. To study the effects of heaving motion on two-phase flow dynamics of marine reactor, Ishida et al. mounted a small scale rig on a suspended platform oscillated in vertical direction by hydraulic device [4]. Murata et al. performed a series of single-phase natural circulation tests in a model reactor with rolling motion to investigate effects of the rolling motion on thermal-hydraulic behavior [5]. Tan et al. also mounted a flow loop on a rolling platform and carried out a series of experiments to study single-phase natural circulation flow, heat transfer and two-phase flow instability [6,7]. Yan and Yu mounted a Passive Residual Heat Removal System (PRHRS) apparatus on rolling equipment to investigate operational characteristics of passive residual heat removal system under rolling motion [8]. The experimental results indicate that loop flow rate oscillates with ship motions and two-phase natural circulation instability may overlap flow oscillation induced by ship motions.

In addition to experiments, a lot of work on theoretical analysis and code modifying have been done to study reactor thermal hydraulics under ocean condition [9,10]. Yan et al. modified the RELAP5 code by adding a module calculating the effect of rolling motion and introducing new flow and heat transfer models to investigate the thermal-hydraulics behavior of nuclear power systems [11]. JAERI developed RETRAN-02/GRAV code, which was a modified version of RETRAN-02 code, to simulate the experimental voyages of Mutsu and to analyze the effects of ship motions on natural circulation of deep sea research reactor DRX [12]. RETRAN-03/MOV code for simulating ship reactor under multi-dimensional motions was also developed from RETRAN-03 in Seoul National University [13].

China plans to build FNPP with proprietary intellectual property rights in the 13th Five-Year Plan. A totally new code PNCMC (Program for Natural Circulation under Motion Conditions) written in C++ language is developed to simulate the natural circulation of FNPP, which is also required to easily add new models, processes and components. One-dimensional two-fluid model in noninertial frame of reference is proposed in this study. Staggered mesh, finite volume method and semi-implicit first order upwind scheme are used to discretize and solve mass, momentum and energy equations. PNCMC code is also verified by single-phase and two-phase natural circulation experimental data.

2. Two-fluid model in noninertial frame of reference

When studying the thermal hydraulic of marine reactor that has an accelerated motion relative to the ground, if taking the ground as the frame of reference, the flow of fluid is a combination of the reactor motion relative to the ground and fluid motion relative to system boundary (pipe wall). However, if we build a noninertial frame which has the same motion as reactor, the flow of fluid in this frame is just the fluid motion relative to the system boundary, which is the main concern in most engineering fields. The boundary condition is changed from the moving condition in inertial frame of reference to stationary boundary condition in noninertial frame of reference; so the mathematical description of flow is simplified.

Thermal dynamic variables and net viscous stress are independent of the frame of reference and the continuity and energy equations are unchanged in a noninertial frame of reference. The primary effect of a noninertial frame is the addition of extra body force term to the momentum equations, and the extra body force term named inertial force arises from the motion of the noninertial frame.

Combining dynamics of mass point in noninertial frame with two-fluid model in inertial frame [14], one-dimensional two-fluid model in noninertial frame of reference is proposed as follows:

Gas continuity equation: (1) $$\begin{array}{ccc}& & {\mathit{\rho }}_g\frac{\partial {\mathit{\alpha }}_g}{\partial t}+{\mathit{\alpha }}_g\frac{\partial {\mathit{\rho }}_g}{\partial t}+\frac{\partial }{\partial z}\left({\mathit{\alpha }}_g{\mathit{\rho }}_g{V}_g\right)\hfill \\ & & =-\frac{H_{\mathrm{ig}}\left(T_s-T_g\right)+H_{\mathrm{if}}\left(T_s-T_f\right)}{h_g^{*}-h_f^{*}}+{\Gamma }_{\mathrm{wg}}\hfill \end{array}$$(1)

Liquid continuity equation: (2) $$\begin{array}{ccc}& & -{\mathit{\rho }}_f\frac{\partial {\mathit{\alpha }}_g}{\partial t}+{\mathit{\alpha }}_f\frac{\partial {\mathit{\rho }}_f}{\partial t}+\frac{\partial }{\partial z}\left({\mathit{\alpha }}_f{\mathit{\rho }}_f{V}_f\right)\hfill \\ & & =\frac{H_{\mathrm{ig}}\left(T_s-T_g\right)+H_{\mathrm{if}}\left(T_s-T_f\right)}{h_g^{*}-h_f^{*}}-{\Gamma }_{\mathrm{wg}}\hfill \end{array}$$(2) where h*_g, h*_f associated with bulk interface mass transfer are defined as(3) $$\begin{array}{ccc}\hfill h_g^{*}& =& 0.5\left[\left(h_{\mathrm{sg}}+h_g\right)+\eta \left(h_{\mathrm{sg}}-h_g\right)\right],\hfill \\ \hfill h_f^{*}& =& 0.5\left[\left(h_{\mathrm{sf}}+h_f\right)-\eta \left(h_{\mathrm{sf}}-h_f\right)\right]\hfill \end{array}$$(3) η = 1 for Γ_ig ⩾ 0 and η = −1 for Γ_ig < 0.

Gas momentum equation: (4) $$\begin{array}{ccc}& & {\alpha }_g{\rho }_g\frac{\partial V_g}{\partial t}+{\alpha }_g{\rho }_g{V}_g\frac{\partial V_g}{\partial z}\hfill \\ & & =-{\alpha }_g\frac{\partial P}{\partial z}+{\Gamma }_g\left(V_i-V_g\right)-p_w{\alpha }_{\mathrm{gw}}{\tau }_{\mathrm{gw}}+M_g^d\hfill \\ & & +⟨{\alpha }_g{\rho }_g\left(g_z+f_{\mathrm{ine},z}\right)⟩\hfill \end{array}$$(4)

Liquid momentum equation: (5) $$\begin{array}{ccc}& & {\alpha }_f{\rho }_f\frac{\partial V_f}{\partial t}+{\alpha }_f{\rho }_f{V}_f\frac{\partial V_f}{\partial z}\hfill \\ & & =-{\alpha }_f\frac{\partial P}{\partial z}+{\Gamma }_f\left(V_i-V_f\right)-p_w{\alpha }_{\mathrm{fw}}{\tau }_{\mathrm{fw}}-M_g^d\hfill \\ & & +⟨{\alpha }_f{\rho }_f\left(g_z+f_{\mathrm{ine},z}\right)⟩\hfill \end{array}$$(5)

Gas energy equation: (6) $$\begin{array}{ccc}& & {\alpha }_g{\rho }_g\frac{\partial U_g}{\partial t}+\left({\rho }_g{U}_g+P\right)\frac{\partial {\alpha }_g}{\partial t}+{\alpha }_g{U}_g\frac{\partial {\rho }_g}{\partial t}\hfill \\ & & +\frac{\partial }{\partial z}\left({\alpha }_g{\rho }_g{U}_g{V}_g\right)+P\frac{\partial }{\partial z}\left({\alpha }_g{V}_g\right)\hfill \\ & & =-\frac{h_f^{*}}{h_g^{*}-h_f^{*}}{H}_{\mathrm{ig}}\left(T_s-T_g\right)-\left(\frac{h_g^{*}}{h_g^{*}-h_f^{*}}\right)\hfill \\ & & \times H_{if}\left(T_s-T_f\right)+{\Gamma }_{\mathrm{wg}}\left(\frac{1+\epsilon }{2}{h}_g^{\text{'}}+\frac{1-\epsilon }{2}{h}_f^{\text{'}}\right)\hfill \\ & & +H_{\mathrm{wg}}\left(T_w-T_g\right)+{\Phi }_g\hfill \end{array}$$(6)

Liquid energy equation: (7) $$\begin{array}{ccc}& & {\alpha }_f{\rho }_f\frac{\partial U_f}{\partial t}+\left({\rho }_f{U}_f+P\right)\frac{\partial {\alpha }_f}{\partial t}+{\alpha }_f{U}_f\frac{\partial {\rho }_f}{\partial t}\hfill \\ & & +\frac{\partial }{\partial z}\left({\alpha }_f{\rho }_f{U}_f{V}_f\right)+P\frac{\partial }{\partial z}\left({\alpha }_f{V}_f\right)\hfill \\ & & =\frac{h_f^{*}}{h_g^{*}-h_f^{*}}{H}_{\mathrm{ig}}\left(T_s-T_g\right)+\left(\frac{h_g^{*}}{h_g^{*}-h_f^{*}}\right)H_{if}\left(T_s-T_f\right)\hfill \\ & & -{\Gamma }_{wg}\left(\frac{1+\epsilon }{2}{h}_g^{\text{'}}+\frac{1-\epsilon }{2}{h}_f^{\text{'}}\right)+H_{\mathrm{fw}}\left(T_w-T_f\right)\hfill \\ & & +{\Phi }_f\hfill \end{array}$$(7) where ξ_h is heating perimeter, Ais the cross area of flow channel, T_w is wall temperature, H_wg, H_wf is volumetric wall heat transfer coefficient for gas and for liquid, respectively, Φ is phase energy dissipation terms. h′_g, h′_f associated with near wall mass transfer are defined as (8) $$\begin{array}{ccc}\hfill h_g^{\text{'}}& =& 0.5\left[\left(h_{\mathrm{sg}}+h_g\right)+\epsilon \left(h_{\mathrm{sg}}-h_g\right)\right],\hfill \\ \hfill h_f^{\text{'}}& =& 0.5\left[\left(h_{\mathrm{sf}}+h_f\right)-\epsilon \left(h_{\mathrm{sf}}-h_f\right)\right]\hfill \end{array}$$(8) ϵ = 1 for Γ_wg ⩾ 0 and ϵ = −1 for Γ_wg < 0.

The constitutive relations include models for defining flow regimes and flow regime-related models for total interfacial shear force M^d_g, wall friction $p_w{\mathit{\alpha }}_{\mathrm{fw}}{\mathit{\tau }}_{\mathrm{fw}}$ and $p_w{\mathit{\alpha }}_{\mathrm{gw}}{\mathit{\tau }}_{\mathrm{gw}}$, wall heat transfer coefficient H_wg and H_wf, interphase heat transfer coefficient H_ig and H_if. In PNCMC, the constitutive relations are in accordance with that used in RELAP5/MOD3.3, and the detailed description of those can be found in RELAP5/MOD3 code manual volume IV [15].

3. Inertial forces arise from ship motion

As shown in Figure 1, the ship motion can be decomposed into swaying, pitching, yawing, heaving, rolling and surging. The axis and amplitude of rolling and pitching can be given by ship design compartment. Among these six modes of motion, pitching, heaving and rolling have significant impact on reactor natural circulation, which are the emphasis of our research. Most often, we tend to study pitching and rolling separately, that is because intense pitching and rolling are not coexisting.

Figure 1. Ship movement – the six degrees of freedom.

The inertial forces of a fluid particle are obtained in this general expression, (9) $$\begin{array}{c}\hfill m·{\mathit{f}}_{\mathrm{ine}}=m\left({\mathit{a}}_0+\stackrel{•}{\mathit{\omega }}\times \mathit{r}+\mathit{\omega }\times \left(\mathit{\omega }\times \mathit{r}\right)+2\mathit{\omega }\times \mathit{V}\right)\end{array}$$(9)

The four terms in the brackets of Equation (9)(9) $$\begin{array}{c}\hfill m·{\mathit{f}}_{\mathrm{ine}}=m\left({\mathit{a}}_0+\stackrel{•}{\mathit{\omega }}\times \mathit{r}+\mathit{\omega }\times \left(\mathit{\omega }\times \mathit{r}\right)+2\mathit{\omega }\times \mathit{V}\right)\end{array}$$(9) are all significant. The first term a₀ accounts for the acceleration of noninertial frame relative to stationary frame. It provides the apparent force that pushes occupants back into their seats or makes them tighten their grip on a handrail when a vehicle accelerates. The second term $\stackrel{•}{\omega }\times \mathbf{r}$ is the acceleration caused by angular acceleration of the noninertial frame, so it is of little importance for geophysical flows or for flows in machinery that rotate at a constant rate about a fixed axis. However, it does play a role when rotation speed or the direction of rotation varies with time. The third term ω × (ω × r) is the centripetal acceleration, which depends strongly on the rotation rate and the distance of the fluid particle from the axis of rotation. The final term 2ω × V is the Coriolis acceleration, depending on the fluid particle's velocity, not on its position.

In Figure 2, one noninertial frame OXYZ with coordinate unit vector i j k and one stationary frame OX_sY_sZ_s is built. The origin O is set on the rolling axis. At initial stage, noninertial frame OXYZ and stationary frame OX_sY_sZ_s are overlapping, the unit vector is straight up and contrary to gravity. Noninertial frame OXYZ has a rolling motion around axis S = (S_x, S_y, S_z) and acceleration a₀ = (a_{0, x}, a_{0, y}, a_{0, z}) relative to stationary frame. Σ is an arbitrary cross section of channel and points in Σ belong to the domain {x₀ ± Δx, y₀ ± Δy, z₀ ± Δz}. The unit vector n = (cos θ, cos β, cos γ) represents for positive main stream of channel in noninertial frame OXYZ.

Figure 2. Geometry showing flow channel and noninertial frame.

The rolling of ship is similar to sine wave and is expressed as (10a) $$\phi ={\phi }_{m}sin\left({\varpi }_{d}t\right),$$(10a) (10b) $$\omega ={\phi }_m{\varpi }_{d}cos\left({\varpi }_{d}t\right),$$(10b) (10c) $$\stackrel{•}{\omega }=-{\phi }_m{\varpi }_d^{2}sin\left({\varpi }_{d}t\right),$$(10c) where φ is rolling angle, φ_m is the rolling amplitude and ϖ_d is the rolling frequency.

According to Equation (9)(9) $$\begin{array}{c}\hfill m·{\mathit{f}}_{\mathrm{ine}}=m\left({\mathit{a}}_0+\stackrel{•}{\mathit{\omega }}\times \mathit{r}+\mathit{\omega }\times \left(\mathit{\omega }\times \mathit{r}\right)+2\mathit{\omega }\times \mathit{V}\right)\end{array}$$(9) , f_ine is obtained as (11) $$\begin{array}{ccc}\hfill {\mathit{f}}_{\mathrm{ine}}& =& [a_{0,x}-{\omega }^2\left(S_x{S}_{y}y-S_y^{2}x-S_z^{2}x+S_x{S}_{z}z\right)\hfill \\ & & -\stackrel{•}{\omega }{S}_{y}z+\stackrel{•}{\omega }{S}_{z}y-2\omega S_y{V}_z+2\omega S_z{V}_y]\mathit{i}\hfill \\ & & +[a_{0,y}-{\omega }^2\left(S_y{S}_{z}z-S_z^{2}y-S_x^{2}y+S_x{S}_{y}x\right)\hfill \\ & & -\stackrel{•}{\omega }{S}_{z}x\stackrel{•}{+\omega }{S}_{x}z-2\omega S_z{V}_x+2\omega S_x{V}_z]\mathit{j}.\hfill \\ & & +[a_{0,z}-{\omega }^2\left(S_x{S}_{z}x-S_x^{2}z-S_y^{2}z+S_y{S}_{z}y\right)\hfill \\ & & -\stackrel{•}{\omega }{S}_{x}y+\stackrel{•}{\omega }{S}_{y}x-2\omega S_x{V}_y+2\omega S_y{V}_x]\mathit{k}\hfill \end{array}$$(11)

The flow field in a channel is actually three-dimensional, but the velocity of main stream is much greater than that of cross section, which is suitable for most flow channels. The Coriolis acceleration is perpendicular to velocity, so the component of Coriolis acceleration along main stream is ignored in PNCMC. Then f_{ine, n} is simplified as (12) $$\begin{array}{ccc}\hfill f_{\mathrm{ine},n}& =& [a_{0,x}-{\omega }^2\left(S_x{S}_{y}y-S_y^{2}x-S_z^{2}x+S_x{S}_{z}z\right)\hfill \\ & & -\stackrel{•}{\omega }{S}_{y}z+\stackrel{•}{\omega }{S}_{z}y]cos\theta +[a_{0,y}-{\omega }^2(S_y{S}_{z}z-S_z^{2}y\hfill \\ & & -S_x^{2}y+S_x{S}_{y}x)-\stackrel{•}{\omega }{S}_{z}x\stackrel{•}{+\omega }{S}_{x}z]cos\beta .\hfill \\ & & +[a_{0,z}-{\omega }^2\left(S_x{S}_{z}x-S_x^{2}z-S_y^{2}z+S_y{S}_{z}y\right)\hfill \\ & & -\stackrel{•}{\omega }{S}_{x}y+\stackrel{•}{\omega }{S}_{y}x]cos\gamma \hfill \end{array}$$(12)

For reactor systems and system experimental facilities, the pipeline cross-sectional scale (Δx,Δy,Δz) is much smaller than system spatial scale (x₀,y₀,z₀); so we can assume that (13) $$\begin{array}{c}\hfill \frac{\Delta x}{x_0}<<1,\frac{\Delta y}{y_0}<<1,\frac{\Delta z}{z_0}<<1.\end{array}$$(13)

So f_{ine, n} is approximately uniform in cross section Σ.

Then ⟨f_{ine, n}⟩ is written as (14) $$\begin{array}{ccc}\hfill ⟨f_{\mathrm{ine},n}⟩& =& [a_{0,x}-{\omega }^2\left(S_x{S}_y{y}_0-S_y^2{x}_0-S_z^2{x}_0+S_x{S}_z{z}_0\right)\hfill \\ & & -\stackrel{•}{\omega }{S}_y{z}_0+\stackrel{•}{\omega }{S}_z{y}_0]cos\theta +[a_{0,y}\hfill \\ & & -{\omega }^2\left(S_y{S}_z{z}_0-S_z^2{y}_0-S_x^2{y}_0+S_x{S}_y{x}_0\right)\hfill \\ & & -\stackrel{•}{\omega }{S}_z{x}_0+\stackrel{•}{\omega }{S}_x{z}_0]cos\beta .+[a_{0,z}\hfill \\ & & -{\omega }^2(S_x{S}_z{x}_0-S_x^2{z}_0-S_y^2{z}_0+S_y{S}_z{y}_0)\hfill \\ & & -\stackrel{•}{\omega }{S}_x{y}_0+\stackrel{•}{\omega }{S}_y{x}_0]cos\gamma \hfill \end{array}$$(14)

⟨f_{ine, n}⟩ is changed to ⟨f_{ine, z}⟩ in Equations (4)(4) $$\begin{array}{ccc}& & {\alpha }_g{\rho }_g\frac{\partial V_g}{\partial t}+{\alpha }_g{\rho }_g{V}_g\frac{\partial V_g}{\partial z}\hfill \\ & & =-{\alpha }_g\frac{\partial P}{\partial z}+{\Gamma }_g\left(V_i-V_g\right)-p_w{\alpha }_{\mathrm{gw}}{\tau }_{\mathrm{gw}}+M_g^d\hfill \\ & & +⟨{\alpha }_g{\rho }_g\left(g_z+f_{\mathrm{ine},z}\right)⟩\hfill \end{array}$$(4) and (5(5) $$\begin{array}{ccc}& & {\alpha }_f{\rho }_f\frac{\partial V_f}{\partial t}+{\alpha }_f{\rho }_f{V}_f\frac{\partial V_f}{\partial z}\hfill \\ & & =-{\alpha }_f\frac{\partial P}{\partial z}+{\Gamma }_f\left(V_i-V_f\right)-p_w{\alpha }_{\mathrm{fw}}{\tau }_{\mathrm{fw}}-M_g^d\hfill \\ & & +⟨{\alpha }_f{\rho }_f\left(g_z+f_{\mathrm{ine},z}\right)⟩\hfill \end{array}$$(5) ).

During rolling, g_z changes with the angle between channel and vertical direction in stationary frame OX_sY_sZ_s. In some sense, rolling is a rotation with varying speed around a fixed axis, the new unit vector of main stream n′ = (cos θ′, cos β′, cos γ′) at the rolling angle φ in stationary frame OX_sY_sZ_s can be obtained through the steps as follows: $$\begin{array}{ccc}\hfill \stackrel{\wedge }{\mathit{A}}& =& \left[\begin{array}{ccc}\hfill S_x{S}_x& \hfill S_x{S}_y& \hfill S_x{S}_z\\ \hfill S_y{S}_x& \hfill S_y{S}_y& \hfill S_y{S}_z\\ \hfill S_z{S}_x& \hfill S_z{S}_y& \hfill S_z{S}_z\end{array}\right],\stackrel{*}{\mathit{A}}=\left[\begin{array}{ccc}0& -S_z& S_y\\ S_z& 0& -S_x\\ -S_y& S_x& 0\end{array}\right],\hfill \\ \hfill \mathit{I}& =& \left[\begin{array}{ccc}\hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 1& \hfill 0\\ \hfill 0& \hfill 0& \hfill 1\end{array}\right]\hfill \end{array}$$(15) $$\mathit{M}=\stackrel{\wedge }{\mathit{A}}+cos\phi ·\left(\mathit{I}-\stackrel{\wedge }{\mathit{A}}\right)+sin\phi ·\stackrel{*}{\mathit{A}}$$(15) (16) $${\mathit{n}}^{\text{'}}=\mathit{n}·{\mathit{M}}^T$$(16) (17) $$g_z=g·cos{\gamma }^{\text{'}}$$(17) where M^T is the transposition of M, and g is the gravity.

4. Numerical solution

PNCMC uses a staggered mesh shown in Figure 3, in which the velocities V_g, V_f are defined at the mass and energy cell interfaces and the pressure p, gas volume fraction α_g, temperature T, internal energy U and density ρ are defined at the mass, energy cell centers. The difference equations for mass and energy cell are obtained by integrating the mass and energy equations with respect to the spatial variable z, from the junction at z_j to z_{j + 1}. The momentum equations are integrated with respect to the spatial variable from cell center z_k to adjoining cell center z_L.

Figure 3. Schematic of staggered mesh.

Completion of the problem definition requires a choice of two independent variables from the four thermodynamic variables: pressure, temperature, density and specific internal energy. Density is not a good choice because a small error in a solution for density can turn into a significant error in pressure. If pressure is chosen as an independent variable, a small error in the solution for pressure results in an even smaller fractional error in density and temperature. The heat transfer at the interfaces of liquid–gas and fluid–heat structures depends on the temperature differences that must be evaluated implicitly in in the numerical equations. So p, T_f, T_g, α_g, V_f, V_g are chosen as independent variables in PNCMC. Under single liquid phase condition (α_g ⩽ 1.0^{− 6}), the gas temperature is set as local saturation temperature; and under single gas phase condition (α_g ⩾ 0.99999), the liquid temperature is also set as local saturation temperature. When phase change occurs (from single-phase to two-phase) in the control volume, the α_g changes first.

Equations of state provide density and internal energy as functions of pressure and temperature. The relationships are generally nonlinear; therefore, Taylor series expansions (18)(18) $${\rho }_k^{n+1}={\rho }_k^n+{\left(\frac{\partial {\rho }_k}{\partial T_k}\right)}^n\frac{T_k^{n+1}-T_k^n}{\Delta t}+{\left(\frac{\partial {\rho }_k}{\partial P}\right)}^n\frac{P^{n+1}-P^n}{\Delta t},$$(18) –(20(20) $$T_s^{n+1}=T_s^n+{\left(\frac{\mathrm{d}{T}_s}{\mathrm{d}P}\right)}^n\frac{P^{n+1}-P^n}{\Delta t},$$(20) ) about old time values are used to obtain a linear formulation in the new time variables. (18) $${\rho }_k^{n+1}={\rho }_k^n+{\left(\frac{\partial {\rho }_k}{\partial T_k}\right)}^n\frac{T_k^{n+1}-T_k^n}{\Delta t}+{\left(\frac{\partial {\rho }_k}{\partial P}\right)}^n\frac{P^{n+1}-P^n}{\Delta t},$$(18) (19) $$U_k^{n+1}=U_k^n+{\left(\frac{\partial U_k}{\partial T_k}\right)}^n\frac{T_k^{n+1}-T_k^n}{\Delta t}+{\left(\frac{\partial U_k}{\partial P}\right)}^n\frac{P^{n+1}-P^n}{\Delta t},$$(19) (20) $$T_s^{n+1}=T_s^n+{\left(\frac{\mathrm{d}{T}_s}{\mathrm{d}P}\right)}^n\frac{P^{n+1}-P^n}{\Delta t},$$(20) where n denotes known values at the current level in time, n + 1 denotes a value at the next level in time, Δt is the increment in time between levels n and n + 1.

To obtain linear difference equations, total interfacial shear force M^d_g, wall friction p_wα_fwτ_fw and p_wα_gwτ_gw are rearranged in new forms: (21a) $$M_g^d=\mathrm{FI}\left(V_g-V_f\right),$$(21a) (21b) $$p_w{\alpha }_{\mathrm{gw}}{\tau }_{\mathrm{gw}}=\mathrm{FWG}·V_g,$$(21b) (21c) $$p_w{\alpha }_{\mathrm{fw}}{\tau }_{\mathrm{fw}}=\mathrm{FWF}·V_f.$$(21c)

The semi-implicit, one order up-wind difference scheme is applied to obtain difference Equations (22)(22) $$\begin{array}{ccc}& & V_L\left({\rho }_g^n\frac{{\alpha }_g^{n+1}-{\alpha }_g^n}{\Delta t}+{\alpha }_g^n\frac{{\rho }_g^{n+1}-{\rho }_g^n}{\Delta t}\right)+{\left({\stackrel{•}{\alpha }}_g{\stackrel{•}{\rho }}_g\right)}_{j+1}^n\hfill \\ & & \times A_{j+1}{V}_{g,j+1}^{n+1}-{\left({\stackrel{•}{\alpha }}_g{\stackrel{•}{\rho }}_g\right)}_j^n{A}_j{V}_g^{n+1}\hfill \\ & & =-V_L\frac{H_{\mathrm{ig}}^n\left(T_s^{n+1}-T_g^{n+1}\right)+H_{if}^n\left(T_s^{n+1}-T_f^{n+1}\right)}{{\left(h_g^{*}-h_f^{*}\right)}^n}\hfill \\ & & +V_L{\Gamma }_{\mathrm{wg}}^n.\hfill \end{array}$$(22) –(27(27) $$\begin{array}{ccc}& & {\left({\alpha }_f{\rho }_f\right)}_j^n\Delta x_j\frac{V_{f,j}^{n+1}-V_{f,j}^n}{\Delta t}+\frac{1}{2}{\left(\stackrel{•}{{\alpha }_f}\stackrel{•}{{\rho }_f}\right)}_j^n[{\left(V_f\right)}_L^n{\left(\stackrel{•}{V_f}\right)}_L^n\hfill \\ & & -{\left(V_f\right)}_K^n{\left(\stackrel{•}{V_f}\right)}_K^n]=-{\left({\alpha }_f\right)}_j^n\left(P_L^{n+1}-P_K^{n+1}\right)\hfill \\ & & +{\left({\Gamma }_f\right)}_j^n\left(V_i^{n+1}-V_{f,j}^{n+1}\right)\Delta x_j-{\left(\stackrel{•}{{\alpha }_f}\stackrel{•}{{\rho }_f}\right)}_j^n\hfill \\ & & \times C_{\mathrm{loss}}{V}_{f,j}^n{V}_{f,j}^{n+1}+{\mathrm{FI}}_j^n\left(V_{g,j}^{n+1}-V_{f,j}^{n+1}\right)\Delta x_j\hfill \\ & & -{\mathrm{FWF}}_j^{n+1}{V}_{f,j}^{n+1}\Delta x_j+{\left({\mathit{\alpha }}_f{\mathit{\rho }}_f\right)}_j^n{\left(g_z+f_{\mathrm{ine},z}\right)}_j^{n+1}\Delta x_j\hfill \end{array}$$(27) ).

The gas continuity difference equation is (22) $$\begin{array}{ccc}& & V_L\left({\rho }_g^n\frac{{\alpha }_g^{n+1}-{\alpha }_g^n}{\Delta t}+{\alpha }_g^n\frac{{\rho }_g^{n+1}-{\rho }_g^n}{\Delta t}\right)+{\left({\stackrel{•}{\alpha }}_g{\stackrel{•}{\rho }}_g\right)}_{j+1}^n\hfill \\ & & \times A_{j+1}{V}_{g,j+1}^{n+1}-{\left({\stackrel{•}{\alpha }}_g{\stackrel{•}{\rho }}_g\right)}_j^n{A}_j{V}_g^{n+1}\hfill \\ & & =-V_L\frac{H_{\mathrm{ig}}^n\left(T_s^{n+1}-T_g^{n+1}\right)+H_{if}^n\left(T_s^{n+1}-T_f^{n+1}\right)}{{\left(h_g^{*}-h_f^{*}\right)}^n}\hfill \\ & & +V_L{\Gamma }_{\mathrm{wg}}^n.\hfill \end{array}$$(22)

The liquid continuity difference equation is (23) $$\begin{array}{ccc}& & V_L\left(-{\rho }_f^n\frac{{\alpha }_g^{n+1}-{\alpha }_g^n}{\Delta t}+{\alpha }_f^n\frac{{\rho }_f^{n+1}-{\rho }_f^n}{\Delta t}\right)+{\left({\stackrel{•}{\alpha }}_f{\stackrel{•}{\rho }}_f\right)}_{j+1}^n\hfill \\ & & \times A_{j+1}{V}_{f,j+1}^{n+1}-{\left({\stackrel{•}{\alpha }}_f{\stackrel{•}{\rho }}_f\right)}_j^n{A}_j{V}_f^{n+1}\hfill \\ & & =+V_L\frac{H_{\mathrm{ig}}^n\left(T_s^{n+1}-T_g^{n+1}\right)+H_{if}^n\left(T_s^{n+1}-T_f^{n+1}\right)}{{\left(h_g^{*}-h_f^{*}\right)}^n}\hfill \\ & & -V_L{\Gamma }_{\mathrm{wg}}^n.\hfill \end{array}$$(23)

The gas energy difference equation is (24) $$\begin{array}{ccc}& & V_L[{\left({\alpha }_g{\rho }_g\right)}^n\frac{U_g^{n+1}-U_g^n}{\Delta t}+{\left({\alpha }_g{U}_g\right)}^n\frac{{\rho }_g^{n+1}-{\rho }_g^n}{\Delta t}\hfill \\ & & +{\left({\alpha }_g{U}_g+P\right)}^n\frac{{\alpha }_g^{n+1}-{\alpha }_g^n}{\Delta t}]+{\left(\stackrel{•}{{\alpha }_g}\stackrel{•}{{\rho }_g}\stackrel{•}{U_g}\right)}_{j+1}^n\hfill \\ & & \times A_{j+1}{V}_{g,j+1}^{n+1}-{\left(\stackrel{•}{{\alpha }_g}\stackrel{•}{{\rho }_g}\stackrel{•}{U_g}\right)}_j^n{A}_j{V}_{g,j}^{n+1}\hfill \\ & & =-P^n\left[{\stackrel{•}{\alpha }}_{g,j+1}^n{A}_{j+1}{V}_{g,j+1}^{n+1}-{\stackrel{•}{\alpha }}_{g,j}^n{A}_j{V}_g^{n+1}\right]\hfill \\ & & -V_L{\left(\frac{h_f^{*}}{h_g^{*}-h_f^{*}}\right)}^n{H}_{\mathrm{ig}}^n\left(T_s^{n+1}-T_g^{n+1}\right)\hfill \\ & & -V_L{\left(\frac{h_g^{*}}{h_g^{*}-h_f^{*}}\right)}^n{H}_{\mathrm{if}}^n\left(T_s^{n+1}-T_f^{n+1}\right)\hfill \\ & & +V_L{\Gamma }_{\mathrm{wg}}^n{\left(\frac{1+\epsilon }{2}{h}_g^{\text{'}}+\frac{1-\epsilon }{2}{h}_f^{\text{'}}\right)}^n\hfill \\ & & +H_{\mathrm{wg}}^n{V}_L(T_w^{n+1}-T_g^{n+1})\hfill \end{array}$$(24)

The liquid energy difference equation is (25) $$\begin{array}{ccc}& & V_L[{\left({\alpha }_f{\rho }_f\right)}^n\frac{U_f^{n+1}-U_f^n}{\Delta t}+{\left({\alpha }_f{U}_f\right)}^n\frac{{\rho }_f^{n+1}-{\rho }_f^n}{\Delta t}\hfill \\ & & -{\left({\alpha }_f{U}_f+P\right)}^n\frac{{\alpha }_g^{n+1}-{\alpha }_g^n}{\Delta t}]+{\left(\stackrel{•}{{\alpha }_f}\stackrel{•}{{\rho }_f}\stackrel{•}{U_f}\right)}_{j+1}^n\hfill \\ & & \times A_{j+1}{V}_{f,j+1}^{n+1}-{\left(\stackrel{•}{{\alpha }_f}\stackrel{•}{{\rho }_f}\stackrel{•}{U_f}\right)}_j^n{A}_j{V}_{f,j}^{n+1}\hfill \\ & & =-P^n\left[{\stackrel{•}{\alpha }}_{f,j+1}^n{A}_{j+1}{V}_{f,j+1}^{n+1}-{\stackrel{•}{\alpha }}_{f,j}^n{A}_j{V}_{f,j}^{n+1}\right]\hfill \\ & & +V_L{\left(\frac{h_f^{*}}{h_g^{*}-h_f^{*}}\right)}^n{H}_{\mathrm{ig}}^n\left(T_s^{n+1}-T_g^{n+1}\right)\hfill \\ & & +V_L{\left(\frac{h_g^{*}}{h_g^{*}-h_f^{*}}\right)}^n{H}_{\mathrm{if}}^n\left(T_s^{n+1}-T_f^{n+1}\right)\hfill \\ & & -V_L{\Gamma }_{\mathrm{wg}}^n{\left(\frac{1+\epsilon }{2}{h}_g^{\text{'}}+\frac{1-\epsilon }{2}{h}_f^{\text{'}}\right)}^n\hfill \\ & & +H_{\mathrm{wf}}^{n+1}{V}_L\left(T_w^{n+1}-T_f^{n+1}\right)\hfill \end{array}$$(25)

The gas momentum difference equation is (26) $$\begin{array}{ccc}& & {\left({\alpha }_g{\rho }_g\right)}_j^n\Delta x_j\frac{V_{g,j}^{n+1}-V_{g,j}^n}{\Delta t}+\frac{1}{2}{\left(\stackrel{•}{{\alpha }_g}\stackrel{•}{{\rho }_g}\right)}_j^n[{\left(V_g\right)}_L^n{\left(\stackrel{•}{V_g}\right)}_L^n\hfill \\ & & -{\left(V_g\right)}_K^n{\left(\stackrel{•}{V_g}\right)}_K^n]=-{\left({\alpha }_g\right)}_j^n\left(P_L^{n+1}-P_K^{n+1}\right)\hfill \\ & & +{\left({\Gamma }_g\right)}_j^n\left(V_i^{n+1}-V_{g,j}^{n+1}\right)\Delta x_j-{\left(\stackrel{•}{{\alpha }_g}\stackrel{•}{{\rho }_g}\right)}_j^n{C}_{\mathrm{loss}}{V}_{g,j}^n{V}_{g,j}^{n+1}\hfill \\ & & -F{I}_j^n\left(V_{g,j}^{n+1}-V_{f,j}^{n+1}\right)\Delta x_j-{\mathrm{FWG}}_j^{n+1}{V}_{g,j}^{n+1}\Delta x_j\hfill \\ & & +{\left({\alpha }_g{\rho }_g\right)}_j^n{\left(g_z+f_{\mathrm{ine},z}\right)}_j^{n+1}\Delta x_j\hfill \end{array}$$(26)

The liquid momentum difference equation is (27) $$\begin{array}{ccc}& & {\left({\alpha }_f{\rho }_f\right)}_j^n\Delta x_j\frac{V_{f,j}^{n+1}-V_{f,j}^n}{\Delta t}+\frac{1}{2}{\left(\stackrel{•}{{\alpha }_f}\stackrel{•}{{\rho }_f}\right)}_j^n[{\left(V_f\right)}_L^n{\left(\stackrel{•}{V_f}\right)}_L^n\hfill \\ & & -{\left(V_f\right)}_K^n{\left(\stackrel{•}{V_f}\right)}_K^n]=-{\left({\alpha }_f\right)}_j^n\left(P_L^{n+1}-P_K^{n+1}\right)\hfill \\ & & +{\left({\Gamma }_f\right)}_j^n\left(V_i^{n+1}-V_{f,j}^{n+1}\right)\Delta x_j-{\left(\stackrel{•}{{\alpha }_f}\stackrel{•}{{\rho }_f}\right)}_j^n\hfill \\ & & \times C_{\mathrm{loss}}{V}_{f,j}^n{V}_{f,j}^{n+1}+{\mathrm{FI}}_j^n\left(V_{g,j}^{n+1}-V_{f,j}^{n+1}\right)\Delta x_j\hfill \\ & & -{\mathrm{FWF}}_j^{n+1}{V}_{f,j}^{n+1}\Delta x_j+{\left({\mathit{\alpha }}_f{\mathit{\rho }}_f\right)}_j^n{\left(g_z+f_{\mathrm{ine},z}\right)}_j^{n+1}\Delta x_j\hfill \end{array}$$(27)

The scalar variable $\underset{k}{\overset{•}{\varphi }}$ at junction j is defined as (28) $$\begin{array}{c}\hfill \underset{k,j}{\overset{•}{\varphi }}=0.5({\varphi }_{k,L}+{\varphi }_{k,K})+0.5\frac{\left|V_{k,j}\right|}{V_{k,j}}\left({\varphi }_{k,L}+{\varphi }_{k,K}\right)\end{array}$$(28)

The solving technique of above equation group is taken from RELAP5 [15]. Through a series of matrix transformation, we can obtain an N × N system of linear equations for the new time pressures. (29) $$\begin{array}{c}\hfill \left[\begin{array}{ccc}{A}_{11}& ...& A_{1N}\\ ⋮& & ⋮\\ A_{N1}& \cdots & A_{NN}\end{array}\right]\left[\begin{array}{c}{P}_1^{n+1}\\ ⋮\\ P_2^{n+1}\end{array}\right]=\left[\begin{array}{c}{B}_1\\ ⋮\\ B_N\end{array}\right]\end{array}$$(29)

Successive Over Relaxation (SOR) method is used to solve Equation (29)(29) $$\begin{array}{c}\hfill \left[\begin{array}{ccc}{A}_{11}& ...& A_{1N}\\ ⋮& & ⋮\\ A_{N1}& \cdots & A_{NN}\end{array}\right]\left[\begin{array}{c}{P}_1^{n+1}\\ ⋮\\ P_2^{n+1}\end{array}\right]=\left[\begin{array}{c}{B}_1\\ ⋮\\ B_N\end{array}\right]\end{array}$$(29) and to obtain the new time pressure for each volume. Then, the pressures are substituted in Equations (26(26) $$\begin{array}{ccc}& & {\left({\alpha }_g{\rho }_g\right)}_j^n\Delta x_j\frac{V_{g,j}^{n+1}-V_{g,j}^n}{\Delta t}+\frac{1}{2}{\left(\stackrel{•}{{\alpha }_g}\stackrel{•}{{\rho }_g}\right)}_j^n[{\left(V_g\right)}_L^n{\left(\stackrel{•}{V_g}\right)}_L^n\hfill \\ & & -{\left(V_g\right)}_K^n{\left(\stackrel{•}{V_g}\right)}_K^n]=-{\left({\alpha }_g\right)}_j^n\left(P_L^{n+1}-P_K^{n+1}\right)\hfill \\ & & +{\left({\Gamma }_g\right)}_j^n\left(V_i^{n+1}-V_{g,j}^{n+1}\right)\Delta x_j-{\left(\stackrel{•}{{\alpha }_g}\stackrel{•}{{\rho }_g}\right)}_j^n{C}_{\mathrm{loss}}{V}_{g,j}^n{V}_{g,j}^{n+1}\hfill \\ & & -F{I}_j^n\left(V_{g,j}^{n+1}-V_{f,j}^{n+1}\right)\Delta x_j-{\mathrm{FWG}}_j^{n+1}{V}_{g,j}^{n+1}\Delta x_j\hfill \\ & & +{\left({\alpha }_g{\rho }_g\right)}_j^n{\left(g_z+f_{\mathrm{ine},z}\right)}_j^{n+1}\Delta x_j\hfill \end{array}$$(26) ) and (27(27) $$\begin{array}{ccc}& & {\left({\alpha }_f{\rho }_f\right)}_j^n\Delta x_j\frac{V_{f,j}^{n+1}-V_{f,j}^n}{\Delta t}+\frac{1}{2}{\left(\stackrel{•}{{\alpha }_f}\stackrel{•}{{\rho }_f}\right)}_j^n[{\left(V_f\right)}_L^n{\left(\stackrel{•}{V_f}\right)}_L^n\hfill \\ & & -{\left(V_f\right)}_K^n{\left(\stackrel{•}{V_f}\right)}_K^n]=-{\left({\alpha }_f\right)}_j^n\left(P_L^{n+1}-P_K^{n+1}\right)\hfill \\ & & +{\left({\Gamma }_f\right)}_j^n\left(V_i^{n+1}-V_{f,j}^{n+1}\right)\Delta x_j-{\left(\stackrel{•}{{\alpha }_f}\stackrel{•}{{\rho }_f}\right)}_j^n\hfill \\ & & \times C_{\mathrm{loss}}{V}_{f,j}^n{V}_{f,j}^{n+1}+{\mathrm{FI}}_j^n\left(V_{g,j}^{n+1}-V_{f,j}^{n+1}\right)\Delta x_j\hfill \\ & & -{\mathrm{FWF}}_j^{n+1}{V}_{f,j}^{n+1}\Delta x_j+{\left({\mathit{\alpha }}_f{\mathit{\rho }}_f\right)}_j^n{\left(g_z+f_{\mathrm{ine},z}\right)}_j^{n+1}\Delta x_j\hfill \end{array}$$(27) ) to obtain the new time velocities. At last, pressures and velocities are substituted in Equations (22(22) $$\begin{array}{ccc}& & V_L\left({\rho }_g^n\frac{{\alpha }_g^{n+1}-{\alpha }_g^n}{\Delta t}+{\alpha }_g^n\frac{{\rho }_g^{n+1}-{\rho }_g^n}{\Delta t}\right)+{\left({\stackrel{•}{\alpha }}_g{\stackrel{•}{\rho }}_g\right)}_{j+1}^n\hfill \\ & & \times A_{j+1}{V}_{g,j+1}^{n+1}-{\left({\stackrel{•}{\alpha }}_g{\stackrel{•}{\rho }}_g\right)}_j^n{A}_j{V}_g^{n+1}\hfill \\ & & =-V_L\frac{H_{\mathrm{ig}}^n\left(T_s^{n+1}-T_g^{n+1}\right)+H_{if}^n\left(T_s^{n+1}-T_f^{n+1}\right)}{{\left(h_g^{*}-h_f^{*}\right)}^n}\hfill \\ & & +V_L{\Gamma }_{\mathrm{wg}}^n.\hfill \end{array}$$(22) )–(25(25) $$\begin{array}{ccc}& & V_L[{\left({\alpha }_f{\rho }_f\right)}^n\frac{U_f^{n+1}-U_f^n}{\Delta t}+{\left({\alpha }_f{U}_f\right)}^n\frac{{\rho }_f^{n+1}-{\rho }_f^n}{\Delta t}\hfill \\ & & -{\left({\alpha }_f{U}_f+P\right)}^n\frac{{\alpha }_g^{n+1}-{\alpha }_g^n}{\Delta t}]+{\left(\stackrel{•}{{\alpha }_f}\stackrel{•}{{\rho }_f}\stackrel{•}{U_f}\right)}_{j+1}^n\hfill \\ & & \times A_{j+1}{V}_{f,j+1}^{n+1}-{\left(\stackrel{•}{{\alpha }_f}\stackrel{•}{{\rho }_f}\stackrel{•}{U_f}\right)}_j^n{A}_j{V}_{f,j}^{n+1}\hfill \\ & & =-P^n\left[{\stackrel{•}{\alpha }}_{f,j+1}^n{A}_{j+1}{V}_{f,j+1}^{n+1}-{\stackrel{•}{\alpha }}_{f,j}^n{A}_j{V}_{f,j}^{n+1}\right]\hfill \\ & & +V_L{\left(\frac{h_f^{*}}{h_g^{*}-h_f^{*}}\right)}^n{H}_{\mathrm{ig}}^n\left(T_s^{n+1}-T_g^{n+1}\right)\hfill \\ & & +V_L{\left(\frac{h_g^{*}}{h_g^{*}-h_f^{*}}\right)}^n{H}_{\mathrm{if}}^n\left(T_s^{n+1}-T_f^{n+1}\right)\hfill \\ & & -V_L{\Gamma }_{\mathrm{wg}}^n{\left(\frac{1+\epsilon }{2}{h}_g^{\text{'}}+\frac{1-\epsilon }{2}{h}_f^{\text{'}}\right)}^n\hfill \\ & & +H_{\mathrm{wf}}^{n+1}{V}_L\left(T_w^{n+1}-T_f^{n+1}\right)\hfill \end{array}$$(25) ) to obtain the new time gas volume fraction and temperatures.

A lumped-parameter solution is available to the user in this basic version to solve the coupling between hydrodynamic volumes and solid structures. And outer iteration is designed to obtain accurate solid temperature. First, we assume a provisional structure temperature T^p_w, and substitute T^p_win Equations (24(24) $$\begin{array}{ccc}& & V_L[{\left({\alpha }_g{\rho }_g\right)}^n\frac{U_g^{n+1}-U_g^n}{\Delta t}+{\left({\alpha }_g{U}_g\right)}^n\frac{{\rho }_g^{n+1}-{\rho }_g^n}{\Delta t}\hfill \\ & & +{\left({\alpha }_g{U}_g+P\right)}^n\frac{{\alpha }_g^{n+1}-{\alpha }_g^n}{\Delta t}]+{\left(\stackrel{•}{{\alpha }_g}\stackrel{•}{{\rho }_g}\stackrel{•}{U_g}\right)}_{j+1}^n\hfill \\ & & \times A_{j+1}{V}_{g,j+1}^{n+1}-{\left(\stackrel{•}{{\alpha }_g}\stackrel{•}{{\rho }_g}\stackrel{•}{U_g}\right)}_j^n{A}_j{V}_{g,j}^{n+1}\hfill \\ & & =-P^n\left[{\stackrel{•}{\alpha }}_{g,j+1}^n{A}_{j+1}{V}_{g,j+1}^{n+1}-{\stackrel{•}{\alpha }}_{g,j}^n{A}_j{V}_g^{n+1}\right]\hfill \\ & & -V_L{\left(\frac{h_f^{*}}{h_g^{*}-h_f^{*}}\right)}^n{H}_{\mathrm{ig}}^n\left(T_s^{n+1}-T_g^{n+1}\right)\hfill \\ & & -V_L{\left(\frac{h_g^{*}}{h_g^{*}-h_f^{*}}\right)}^n{H}_{\mathrm{if}}^n\left(T_s^{n+1}-T_f^{n+1}\right)\hfill \\ & & +V_L{\Gamma }_{\mathrm{wg}}^n{\left(\frac{1+\epsilon }{2}{h}_g^{\text{'}}+\frac{1-\epsilon }{2}{h}_f^{\text{'}}\right)}^n\hfill \\ & & +H_{\mathrm{wg}}^n{V}_L(T_w^{n+1}-T_g^{n+1})\hfill \end{array}$$(24) ) and (25(25) $$\begin{array}{ccc}& & V_L[{\left({\alpha }_f{\rho }_f\right)}^n\frac{U_f^{n+1}-U_f^n}{\Delta t}+{\left({\alpha }_f{U}_f\right)}^n\frac{{\rho }_f^{n+1}-{\rho }_f^n}{\Delta t}\hfill \\ & & -{\left({\alpha }_f{U}_f+P\right)}^n\frac{{\alpha }_g^{n+1}-{\alpha }_g^n}{\Delta t}]+{\left(\stackrel{•}{{\alpha }_f}\stackrel{•}{{\rho }_f}\stackrel{•}{U_f}\right)}_{j+1}^n\hfill \\ & & \times A_{j+1}{V}_{f,j+1}^{n+1}-{\left(\stackrel{•}{{\alpha }_f}\stackrel{•}{{\rho }_f}\stackrel{•}{U_f}\right)}_j^n{A}_j{V}_{f,j}^{n+1}\hfill \\ & & =-P^n\left[{\stackrel{•}{\alpha }}_{f,j+1}^n{A}_{j+1}{V}_{f,j+1}^{n+1}-{\stackrel{•}{\alpha }}_{f,j}^n{A}_j{V}_{f,j}^{n+1}\right]\hfill \\ & & +V_L{\left(\frac{h_f^{*}}{h_g^{*}-h_f^{*}}\right)}^n{H}_{\mathrm{ig}}^n\left(T_s^{n+1}-T_g^{n+1}\right)\hfill \\ & & +V_L{\left(\frac{h_g^{*}}{h_g^{*}-h_f^{*}}\right)}^n{H}_{\mathrm{if}}^n\left(T_s^{n+1}-T_f^{n+1}\right)\hfill \\ & & -V_L{\Gamma }_{\mathrm{wg}}^n{\left(\frac{1+\epsilon }{2}{h}_g^{\text{'}}+\frac{1-\epsilon }{2}{h}_f^{\text{'}}\right)}^n\hfill \\ & & +H_{\mathrm{wf}}^{n+1}{V}_L\left(T_w^{n+1}-T_f^{n+1}\right)\hfill \end{array}$$(25) ). After solving the difference Equations (22)(22) $$\begin{array}{ccc}& & V_L\left({\rho }_g^n\frac{{\alpha }_g^{n+1}-{\alpha }_g^n}{\Delta t}+{\alpha }_g^n\frac{{\rho }_g^{n+1}-{\rho }_g^n}{\Delta t}\right)+{\left({\stackrel{•}{\alpha }}_g{\stackrel{•}{\rho }}_g\right)}_{j+1}^n\hfill \\ & & \times A_{j+1}{V}_{g,j+1}^{n+1}-{\left({\stackrel{•}{\alpha }}_g{\stackrel{•}{\rho }}_g\right)}_j^n{A}_j{V}_g^{n+1}\hfill \\ & & =-V_L\frac{H_{\mathrm{ig}}^n\left(T_s^{n+1}-T_g^{n+1}\right)+H_{if}^n\left(T_s^{n+1}-T_f^{n+1}\right)}{{\left(h_g^{*}-h_f^{*}\right)}^n}\hfill \\ & & +V_L{\Gamma }_{\mathrm{wg}}^n.\hfill \end{array}$$(22) –(27(27) $$\begin{array}{ccc}& & {\left({\alpha }_f{\rho }_f\right)}_j^n\Delta x_j\frac{V_{f,j}^{n+1}-V_{f,j}^n}{\Delta t}+\frac{1}{2}{\left(\stackrel{•}{{\alpha }_f}\stackrel{•}{{\rho }_f}\right)}_j^n[{\left(V_f\right)}_L^n{\left(\stackrel{•}{V_f}\right)}_L^n\hfill \\ & & -{\left(V_f\right)}_K^n{\left(\stackrel{•}{V_f}\right)}_K^n]=-{\left({\alpha }_f\right)}_j^n\left(P_L^{n+1}-P_K^{n+1}\right)\hfill \\ & & +{\left({\Gamma }_f\right)}_j^n\left(V_i^{n+1}-V_{f,j}^{n+1}\right)\Delta x_j-{\left(\stackrel{•}{{\alpha }_f}\stackrel{•}{{\rho }_f}\right)}_j^n\hfill \\ & & \times C_{\mathrm{loss}}{V}_{f,j}^n{V}_{f,j}^{n+1}+{\mathrm{FI}}_j^n\left(V_{g,j}^{n+1}-V_{f,j}^{n+1}\right)\Delta x_j\hfill \\ & & -{\mathrm{FWF}}_j^{n+1}{V}_{f,j}^{n+1}\Delta x_j+{\left({\mathit{\alpha }}_f{\mathit{\rho }}_f\right)}_j^n{\left(g_z+f_{\mathrm{ine},z}\right)}_j^{n+1}\Delta x_j\hfill \end{array}$$(27) ), we substitute T^{n + 1}_g, T^{n + 1}_f in Equation (30)(30) $$\begin{array}{ccc}& & V_s{\left({\mathit{\rho }}_s{C}_{p,s}\right)}^n\frac{T_w^{n+1}-T_w^n}{\Delta t}\hfill \\ & & =V_s{\left(q^{{}^{\text{'}\text{'}\text{'}}}\right)}^{n+1}-\sum _i^{i=N}\left\{V_i\right[H_{\mathrm{wf}}^n\left(T_w^{n+1}-T_f^{n+1}\right)\hfill \\ & & +H_{\mathrm{wg}}^n\left(T_w^{n+1}-T_g^{n+1}\right)\left]\right\},\hfill \end{array}$$(30) to obtain T^{n + 1}_w. (30) $$\begin{array}{ccc}& & V_s{\left({\mathit{\rho }}_s{C}_{p,s}\right)}^n\frac{T_w^{n+1}-T_w^n}{\Delta t}\hfill \\ & & =V_s{\left(q^{{}^{\text{'}\text{'}\text{'}}}\right)}^{n+1}-\sum _i^{i=N}\left\{V_i\right[H_{\mathrm{wf}}^n\left(T_w^{n+1}-T_f^{n+1}\right)\hfill \\ & & +H_{\mathrm{wg}}^n\left(T_w^{n+1}-T_g^{n+1}\right)\left]\right\},\hfill \end{array}$$(30) where N is hydrodynamic volume number connected to structure, and V_i is the volume of hydrodynamic volume.

If the results do not satisfy the convergence criterion |T^{n + 1}_w − T_w^p| ⩽ 0.01, then we make T^p_w = T_w^{n + 1} and repeat above processes. A schematic representation of the algorithm is depicted in Figure 4.

Figure 4. Algorithm to solve the two-phase flow equations.

5. Code validation

Code validation should cover a broad range of conditions and phenomena. A lot of verifications of PNCMC for stationary problem, such as subcool boiling and Type-I instability in Chinese nuclear heating reactor, have been done in [16]. Since only few natural circulation experiments under motion conditions were performed in previous studies, we don't have enough experiment data for validation. In this study, single-phase natural circulation experiments under rolling conditions performed in INET (Institute of nuclear and new energy technology of Tsinghua University) and two-phase natural circulation experiments under rolling condition performed by Tan et al. are used for the validation of PNCMC.

5.1. Single-phase natural circulation under rolling conditions

To understand the natural circulation behavior of the integral type reactor under inclination and rolling conditions, a scaled test facility is designed in INET. The test facility is a reactor simulator (natural circulation test loop) mounted on rolling equipment, and shown schematically in Figure 5. The symmetrically arranged test loop consists of three heating channels (from left to right, HC1, HC2 and HC3), two heat exchangers (HEX_A and HEX_B), a riser section, a nitrogen pressurizer, hot legs (Hot leg A and Hot leg B), lower plenum and connecting pipes. The general configuration of the test loop was determined with reference to NHR-200. The natural circulation flow rate was measured at the outlets of each heating channel and downcomer by ultrasonic flow meters, which have the accuracy of 0.5% over the measured flow rate value. Temperature is measured by thermocouples with accuracy of 1% over the temperature range. The rolling begins from the time that natural circulation reaches steady state under vertical condition.

Figure 5. Schematic of the test facility.

Figure 6 shows the results of one case of natural circulation experiment under rolling condition. As was observed, the heat and flow symmetry between heating channels and hot legs was broken. The mass flow rates of heating channels began to oscillate when the rolling started, the oscillation periods of left and right heating channels were the same as those in rolling motion; however, the oscillation period of middle heating channel was half of the period in rolling motion. The same oscillations were also observed in mass flow rates of hot legs with rolling period; however, oscillation amplitudes were much larger than that of heating channels. The phenomenon was similar to that found by Murata et al. [5].

Figure 6. The result of natural circulation experiment under rolling condition.

The node model for test facility is shown in Figure 7, and the distance from rolling axis (red point in Figure 7) to bottom of test loop is 3.2 m. The coordinates of hydrodynamic volumes’ center and unit vector of main flow stream in the noninertial frame are added when establishing the input card. Boundary conditions of the analysis including the inlet water temperature and the outlet pressure of second side of heat exchangers, and the pressure of simulator were constants. The initial condition of transient simulation corresponds to vertical steady-state. The mass flow rate of second side was kept constant during calculation, which was consistent with experiment.

Figure 7. Nodalization of test facility for PNCMC.

The comparisons between experimental results and PNCMC prediction are shown in Figures 8 and 9. The comparisons indicated that PNCMC can figure out the variation trend of mass flow rate. The pressure drop coefficients (heating channels and heat exchangers; places of abrupt area change) used in the input card were measured in different power level in vertical steady-state experiments; so there was a small numerical deviation between experiment and PNCMC.

Figure 8. Comparison between experiment results and PNCMC prediction, experiment case: φ_m = 22.5°, τ = 13 s.

Figure 9 . Comparison between experiment results and PNCMC prediction, experiment case: ${\phi }_m={15}^{\circ },\mathit{\tau }=8\mathrm{s}$.

5.2. Two-phase natural circulation under rolling conditions

Tan et al. experimentally studied the two-phase flow instability of natural circulation under rolling motion condition [7]. The test facility is shown in Figure 10, and is composed of a simple natural circulation loop and a rolling plate. According to the configuration of the test facility, the natural circulation loop was modeled with 56 volumes, as shown in Figure 11. The detailed local pressure drop coefficients, the length and diameter of pipes, and experiments procedure are not available in [6], which are major influence factors of two-phase natural circulation. Through the incomplete information, it is difficult to simulate accurately of Tan's experiments. So, this validation focused on the repetition of phenomenon, not on numerical accuracy.

Figure 10. Schematic of test facility (Figure 2 in Tan et al. [7]).

Figure 11. The nodalization of Tan's test facility for PNCMC.

Figure 12 shows the mass flow rate variation under rolling motion with increasing heating power from 3 to 13 kW. The rolling motion began at 2500 s, and then the mass flow rate oscillated periodically. As the heating power increased, the oscillation amplitude decreased slowly and then enlarged suddenly when the power reached 12.0 kW. Figure 13 indicates that vapor bubbles generated around the trough point of the flow oscillation, and single-phase flow and two-phase flow alternatively appeared in heating pipe, which was also observed in Tan's experiments. This oscillation was similar to density wave oscillation (DWO), but it was not self-excited; the flow would become steady single-phase flow if rolling motion stops. Increasing heating power continuously, the variation of mass flow rate was no longer approximate sine wave. In one rolling period, the mass flow rate increased faster than decreased, and a short-lived increase occurred in the process of decrease marked by arrow in Figure 14. This irregular complex flow oscillation was a transition process, and the oscillation kept evolving until an overlapping of flow oscillation caused by rolling motion and DWO occurred in Figure 15. The prediction of PNCMC was consistent with experimental result. In one rolling period, there were three short periodic oscillations from 3.1 to 3.6 s, which was because the period of DWO was about 3 s. As shown in Figure 11, the overlapping effect was obvious and enhanced the oscillation, and the system became more unstable. To sum up, PNCMC has the ability to simulate two-phase flow under rolling motion.

Figure 12. Mass flow rate variation under rolling motion.

Figure 13. Mass flow rate variation at heating power equals to 6.0 kW.

Figure 14. Irregular complex flow oscillation.

Figure 15. Overlapping of flow oscillation caused by rolling motion and density wave oscillation.

6. Conclusions

In the present work, one-dimensional two-fluid model in noninertial frame of reference is proposed. The body force term in the momentum equation is modified as time dependent function and inertial force induced by ship's three-dimensional motion, such as rolling, heaving and inclination are added into this term. The constitutive relations including models for defining flow regimes and flow regime-related models for total interfacial shear force, wall friction, wall heat transfer and interphase heat transfer are in accordance with that used in RELAP5/MOD3.3. Staggered mesh, finite volume method, semi-implicit first order upwind scheme, lumped-parameter solution of heat structures and SOR method are used to discretize and solve two-phase mass, momentum and energy equations. PNCMC code is developed based on these difference equations and compiled using C++ language.

Single-phase natural circulation experiments under rolling conditions performed in INET and two-phase natural circulation experiments under rolling condition performed by Tan et al. are used for the validation of PNCMC. The calculation results are in good agreement with single-phase experiments. PNCMC can figure out the phenomenon that was found in two-phase experiments. However, much more two-phase natural circulation experiments are needed to evaluate the accuracy of PNCMC.

Nomenclature
A	=	Area of cross section (m2)
a0	=	Acceleration of noninertial frame relative to stationary frame (m/s2)
Cp	=	Specific heat at constant pressure (J/kg · K)
FI	=	Interphase drag coefficients (liquid, vapor) (s−1)
FWF, FWG	=	Wall drag coefficients (liquid, vapor) (s−1)
fine	=	Extra acceleration (m/s2)
g	=	Gravity acceleration (−9.8 m/s2)
Hif	=	Volumetric heat transfer coefficient from interface to liquid (W/K · m3)
Hig	=	Volumetric heat transfer coefficient from interface to gas (W/K · m3)
HLoss	=	Pressure drop (Pa)
h	=	Enthalpy (J/kg)
Hwf	=	volumetric wall heat transfer coefficient for liquid (W/Km3)
Hwg	=	volumetric wall heat transfer coefficient for gas (W/Km3)
Mdg	=	Total interfacial force (N)
m	=	Mass (kg)
P	=	Pressure (Pa)
pw	=	Channel perimeter (m)
q	=	Heating power (w)
q'''	=	Power source (w/m3)
r	=	Radial vector of fluid particle in noninertial frame (m)
T	=	Temperature (K)
t	=	Time (s)
U	=	Specific internal energy (J/kg)
V	=	Velocity in noninertial frame (m/s), volume (m3)
x, y, z	=	Spatial coordinate (m)
Greek letters	=
α	=	Gas volume fraction
αgw	=	Gas volume fraction near wall
αfw	=	liquid volume fraction near wall
ζh	=	heating perimeter (m)
ρ	=	Density (kg/m3)
Γwg	=	Volumetric mass exchange rate near wall (kg/m3 · s)
φ	=	Rolling angle (rad)
ω	=	Angle velocity of noninertial frame (rad/s)
τ	=	Shear stresses (N), rolling period (s)
$\stackrel{•}{\omega }$	=	Angle acceleration of noninertial frame (rad/s2)
ϖd	=	Rolling frequency of noninertial frame (Hz)
Subscripts	=
f	=	Liquid phase
g	=	Gas phase
j	=	Spatial noding indices for junctions
K	=	Spatial noding index for volumes
L	=	Spatial noding index for volumes
k	=	k = f or k = g
m	=	Maximum
i	=	Two-phase interface
s	=	Saturation state, heat structure
sf	=	Saturation state of liquid
sg	=	Saturation state of gas
w	=	Pipe wall
z	=	Main stream
Superscripts	=
n, n + 1	=	Time level index
•	=	Donored quantity
*	=	Bulk/saturation property

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the National Science Fund for Distinguished Young Scholars [No.11325526].