The flow in the reactor satisfies three conservation equations of fluid, namely mass continuity, momentum conservation and energy conservation.

The mass conservation is expressed as ∂ ρ ∂ t + ∇ ⋅ ( ρ v ⇀ ) = 0 (1) The equation of energy conservation is ∂ ∂ t ( ρ E ) + ∇ ⋅ [ v ⇀ ( ρ E + p ) ] = ∇ ⋅ ( k e f f ∇ T + τ ¯ ¯ e f f ⋅ v ⇀ ) (2) The equation of momentum conservation can be expressed as ∂ ∂ t ( ρ v ⇀ ) + ∇ ⋅ ( ρ v ⇀ v ⇀ ) = − ∇ p + ∇ ⋅ [ μ ( ∇ v ⇀ + ∇ v ⇀ T ) ] + ρ g ⇀ + S i (3) where S _i is the momentum source term introduced by the porous medium model, which can be denoted as S i = − ( μ α v i + C 2 1 2 ρ | v | v j ) (4) where 1/α and C₂ are viscous resistance coefficient and inertia resistance coefficient, respectively, which need to be determined by fine mesh modeling of local structure.

Besides, when the reactor is subjected to ocean motions which can produce additional accelerations, those accelerations must be introduced into field equations, as discussed in Oscillating Motion.

Shear Stress Transport (SST) turbulence model is one of the k − ω models, which is proposed and developed by Menter (1994), Menter and Esch (2001), Menter et al. (2003). It uses a blending function to realize a smooth transition between the standard k − ω model near a wall and the k − ω model in the free stream. The basic turbulent shear stress is considered in the SST model to correct the turbulent viscosity. Meanwhile, in order to ensure good prediction in the near-wall region and the far field, the cross diffusion term and the mixing function in the transport equation are also modified. The transport equation is as follows: ρ D k D t = ∂ ∂ x i ( Γ k ∂ k ∂ x i ) + G ˜ k − Y k + S k (5) ρ D ω D t = ∂ ∂ x i ( Γ ω ∂ ω ∂ x i ) + G ω − Y ω + D ω + S ω (6) where k and ω are turbulent kinetic energy and specific turbulent dissipation rate (or turbulent characteristic frequency); Γ k and Γ ω represent diffusion coefficients, and are expressed as Γ k = μ + μ t σ k , Γ ω = μ + μ t σ ω , respectively. The Prandtl number σ and Turbulent viscosity μ are formulated as follows: σ k = 1 F 1 / σ k , 1 + ( 1 − F 1 ) / σ k , 2 (7) σ ω = 1 F 1 / σ ω , 1 + ( 1 − F 1 ) / σ ω , 2 (8) μ t = ρ k ω 1 max [ 1 α , S F 2 α 1 ω ] (9) where the blending fuction F ₁ is defined by F 1 = tanh ( Φ 1 4 ) , with Φ 1 = min [ max ( k 0.09 ω y , 500 μ ρ y 2 ω ) , 4 ρ k σ ω , 2 D ω + y 2 ] , D ω + = max ( 2 ρ 1 σ ω , 2 1 ω ∂ k ∂ x j ∂ ω ∂ x j , 10 − 2 )

And the second blending function F ₂ is described as F 2 = tanh ( Φ 2 2 ) , with Φ 2 = max ( k 0.09 ω y , 500 μ ρ y 2 ω ) .

The turbulence generation terms G ˜ k and G ω are expressed as G ˜ k = min ( G k , 10 ρ β ∗ k ω ) and G ω = α ∞ v t G k , where α ∞ = F 1 / α ∞ , 1 + ( 1 − F 1 ) / α ∞ , 2 , α ∞ , 1 = β i , 1 β ∗ − k 2 σ ω , 1 β ∗ , α ∞ , 2 = β i , 2 β ∗ − k 2 σ ω , 2 β ∗ .

The turbulent dissipation terms Y k and Y ω are expressed as follows Y k = ρ β ∗ k ω , Y ω = ρ β ω 2 . Any constant φ ( σ k ,   σ ω ,   β ,   δ ) of the k − ω SST turbulence model can be obtained from the constant φ 1 ( σ k 1 ,   σ ω 1 ,   β 1 ,   δ 1 ) in the k − ω turbulence model (Wilcox, 1988) and the constant φ 2 ( σ k 2 ,   σ ω 2 ,   β 2 ,   δ 2 ) in the transformed k − ω model (Mohammadi and Pironneau, 1993). The constants of the k − ω SST turbulence model in this study are: σ k , 1 = 1.176 , σ k , 2 = 1.0 , σ ω , 1 = 2.0 , σ ω , 2 = 2.0 , α 1 = 0.31 , β i , 1 = 0.075 , β i , 2 = 0.0828 , and β ∗ = 0.09 .

In this study, porous media model is employed to help simulate the pressure drop and the flow resistance inside the reactor. And the areas where the porous media model is adopted mainly include the core fuel assembly area and the steam generator area.

For single-phase and multiphase media, porous media model can take the form of apparent velocity or physical velocity. Assuming that the fluid is single phase, the velocity of the fluid flowing through the area is called the apparent velocity, which can be described as the flow rate through unit cross-sectional area. The porous media model based on apparent velocity can better simulate the pressure loss in porous media. The resistance of fluid in porous materials is calculated by adding a source term to the Navier-Stokes equation, which determines the velocity and pressure field within a fluid domain. The source term mainly consists of inertia loss term and Darcy viscous resistance term, which can be expressed as: S i = − ( ∑ j = 1 3 D i j μ v j + ∑ j = 1 3 C i j 1 2 ρ | v | v j ) (10) here, C _ij and D _ij are the coefficient matrices of inertia loss and viscous loss, respectively. The pressure drop in the porous media unit is produced by this negative momentum source term.

For simple homogeneous porous media, the diagonal terms 1/α and C2 are introduced into the coefficient matrix C _ij and D _ij , respectively, while the other terms are zero, so S _i is given as follows: S i = − ( μ α v i + C 2 1 2 ρ | v | v j ) (11) where 1/α and C ₂ are viscous resistance coefficient and inertia resistance coefficient, respectively.

When the flow in the porous medium is laminar, the pressure drop is proportional to the velocity, which means the porous medium model can be simplified as Darcy’s Law: Δ p = − μ α v ⇀ (12) Thus, the pressure drop in the three coordinate directions is given as follows: Δ p i = ∑ j = 1 3 μ α i j v j Δ n i (13) where Δn _i represents the thickness of the porous medium in three coordinate directions.

When the speed is relatively high, or when simulating porous plates and tube rows, the viscous loss term sometimes can be ignored, and only the inertial loss term is retained, hence the porous media equation can be simplified to: Δ p = − ∑ j = 1 3 C i j 1 2 ρ | v | v j (14) or given as the pressure drop in three coordinate directions: Δ p i = ∑ j = 1 3 C i j Δ n i 1 2 ρ | v | v j (15) In a porous medium, when the permeability of the medium is large and the geometric scale of the medium does not interact with the scale of the turbulent vortex, it can be considered that the solid matrix has no effect on the generation and dissipation rate of turbulence.

With regard to resistance coefficient, it is generally defined based on the apparent velocity of fluid in porous media. In this study, the inertia resistance coefficient and viscous resistance coefficient are determined by the relationship between pressure drop Δp and velocity v in porous media. Suppose the second-order polynomial fitting relationship between pressure drop and velocity is: Δ p = a 1 v + a 2 v 2 (16) where a ₁ and a ₂ are fitting coefficients. The source term in the N-S equation is the pressure drop per unit length, which is given as follows: S i = − Δ p ⋅ Δ n (17) Then, the resistance coefficients can be obtained bycalculation, which are expressed as follows: C 2 = a 2 ρ Δ n / 2 (18) 1 α = a 1 μ Δ n (19)

In this simulation, the hybrid mesh method is adopted to mesh the established geometric model. As shown in Figures 1A,B, the area of the lower head and the upper chamber is divided by poly-hexa meshes. By contrast, full hexahedral meshes are used in the inlet section, the descending section of the downcommer, the core and the steam generator area, as shown in Figures 1C,D.

Then, the mesh sensitivity is demonstrated by generating meshes of different scales. Under the same cold state condition, the comparison of pressure drop throughout the reactor with different mesh scales is obtained and the desirable mesh scheme for subsequent calculation can be determined.

Table 1 shows the number of meshes based on different meshing schemes. ICEM is employed to establish hexahedral grid for annulus, core and steam generator area, while Fluent Meshing is used for upper and lower head regions for the Poly-hexa Meshing. In order to get a reasonable computing mesh model, different number of meshes is generated in different regions. Then the whole model is verified under standard conditions to determine the computational mesh model needed for subsequent simulation.

Figure 2 compares the pressure drop in different regions of the reactor under different mesh schemes. It can be obtained that the average pressure of each area tends to be stable with the increase of grid number. The results of mesh scheme 3, 4, and 5 show little difference, so mesh scheme 3 is employed in subsequent calculation and analysis.

As discussed before, the flow in the reactor satisfies three conservation equations of fluid, including mass continuity, momentum conservation and energy conservation. When the reactor system is subjected to additional motions, additional forces or accelerations are generated and must be taken into consideration in governing equations.

In this study, considering the effect of oscillating motion, a volume force is introduced into the momentum equation, which is expressed as: ρ ( ∂ u ∂ t + u ⋅ ∇ u ) = − ∇ P + σ + ρ ( g − a ) (20) where g is the gravitational acceleration. The oscillatin axis is 3 m above the bottom of the pressure vessel. During the oscillating process, the decomposition of the gravitational acceleration is required. The additional acceleration a is caused by the oscillating motion, which includes centrifugal force, tangential force, and Coriolis force, and can be described as: a = r × ω ˙ + 2 ω × V + ω × ( ω × r ) (21) where ω ˙ ,   ω ,   V   and   r represent the angular accerelation, angular velocity, fluid velocity and position vector, respectively.

The oscillating angle in the oscillating motion can be derived as (Gao, 1999): φ = φ m ⋅ sin ( ω d t ) (22) where φ m is the maximum oscillating angle and ω d is the average angular velocity, which can be expressed as ω d = 2 π T . Here, T represents the period of oscillating motion.

And the angular velocity and accerelation can be derived as: ω = ω d ⋅ φ m ⋅ cos ( ω d t )   (23) ω ˙ = − ω d ⋅ ω d ⋅ φ m ⋅ sin ( ω d t ) (24) By decomposing the additional acceleration, the expression of the additional motion inertia force is obtained as follows: { F x = − ρ ω 2 ( S x S y y − S y 2 x − S z 2 x + S x S z z ) − ρ ω ˙ S y z + ρ ω ˙ S z y − 2 ρ ω S y W + 2 ρ ω S z V F y = − ρ ω 2 ( S y S z z − S z 2 y − S x 2 y + S x S y x ) − ρ ω ˙ S z x + ρ ω ˙ S x z − 2 ρ ω S z U + 2 ρ ω S x W F z = − ρ ω 2 ( S x S z x − S x 2 z − S y 2 z + S y S z y ) − ρ ω ˙ S x y + ρ ω ˙ S y x − 2 ρ ω S x V + 2 ρ ω S y U (25)

Three dimensional simulation on the cold state flow field is conducted based on a full-scale marine reactor, considering the influence of gravity and oscillating inertia force. And the transient hydrodynamic analysis of the whole model is carried out. The temperature of water used as coolant is set at 230°C. The energy exchange between the core and steam generator is neglected.

In this research, porous media model is adopted in the core fuel assembly area and the steam generator area. In order to simulate the flow field, the resistance coefficient need to be determined according to the relationship between pressure drop and velocity.

The flow distribution plate and the lower head will directly affect the flow distribution into the core and further influence the local flow field, so it is necessary to perform refined solid modeling with the structure details considered. Through the standard model, the distribution of flow rate, velocity as well as pressure is investigated. Then, after introducing additional oscillation, the effect of oscillating angle and elevation on the flow field is studied.

The oscillating models are listed in Table 2 with different oscillating angles and elevations. The completed calculation conditions are consistent with the technical requirements. After the calculation, the flow field and pressure of the core, the lower head, the downcomer and the primary side of the steam generator are presented in the form of streamline and velocity contour. The aim is to explore the influence of oscillating angle and oscillating elevation on the cold hydrodynamics.

The resistance coefficient can be determined by the relationship between pressure drop and velocity in porous media as described previously in Porous Media Model. In the research, the pressure loss of fuel assembly with lattice is calculated without considering the effect of gravity. According to the pressure at the different selected cross-sections, the pressure loss and the flow resistance of the light rod and the lattice area can be obtained.

As shown in Figure 3, the Z1/Z2 cross-section is extracted to calculate the regional resistance coefficient of the lattice, and the pressure difference of Z3/Z4 cross-section is used to calculate the regional resistance coefficient of the light rod.

Both the reactors with and without additional movements are studied. For the standard model, no additional motion is imposed to the reactor in order to show the original flow field. The boundary conditions for the two types of model are explained in the following.

For comparison, the flow field in a static reactor is studied first in Standard Model. No additional motion is imposed to the reactor. The inlet flow distribution, velocity field, and pressure are discussed in this part.

This part focuses on the flow characteristics under oscillating motion. The influene of the oscillating motion is presented by comparing the fluid flow under four different oscillation models, which are introduced in Simulation Model and Conditions. The results of flow field distribution, different oscillating angles and different elevations are analyzed as follows.