Application of an Integral Turbulence Model to Close the Model of an Anisotropic Porous Body as Applied to Rod Structures

Abstract

In practice, often devices are ordered rod structures consisting of a large number of rods. Heat exchangers, fuel assemblies of nuclear reactors, and their cores in the case of using caseless assemblies are examples of such devices. Simulation of heat and mass transfer processes in such devices in porous-body approximation can significantly reduce the required resources compared to computational fluid dynamics (CFD) approaches. The paper describes an integral turbulence model developed for defining anisotropic model parameters of a porous body. The parameters of the integral turbulence model were determined by numerical simulations for assemblies of smooth rods, assemblies with spacer grids, and wire-wrapped fuel assemblies. The results of modeling the flow of a liquid metal coolant in an experimental fuel assembly with local blocking of its flow section in anisotropic porous-body approximation using an integral turbulence model are described. The possibility of using the model of an anisotropic porous body with the integral model of turbulence to describe thermal-hydraulic processes during fluid flow in rod structures is confirmed.

1. Introduction

In engineering practice, integral and local approaches are used to simulate thermal-hydraulic processes in the cores of nuclear power plants and other heat-exchange equipment. The integral approach includes sub-channel and porous body models. These models operate with process averaged characteristics. Local models operate with local characteristics and allow one to determine their value at any point in the computational domain. They are implemented in the so-called computational fluid dynamics (CFD) codes. Both approaches are based on the Navier–Stokes (NS) equations, which describe mass and momentum conservation. In addition to those equations, the energy conservation equation is solved.

Sub-channel models operate with averaged characteristics, which are of the main interest in most problems but allow one to determine their value only at the corresponding nodes of the elements into which the computational domain is divided. These methods were developed [1,2] and continue to be used now. A large number of codes have been developed using this approach. The ability to integrate sub-channel and system code in one software product [3], as well as relatively low requirements for computer resources [4], are the advantages of this approach. The main disadvantage of this approach is the need to determine the closing relations, in particular, the coefficients of interchannel mixing ratio, which directly depend on the corresponding geometry and flow regime [3,4,5,6]. These relations are not directly derived from the three-dimensional equations for the flow and should be determined mainly experimentally for similar conditions.

The porous-body model has lower computational costs compared to CFD codes and is comparable to sub-channel models. In contrast to sub-channel models, the closing relations are obtained by averaging three-dimensional equations [7] and have a clear physical meaning; therefore, they are easier to determine, and they work in a wide range of geometric and operating parameters. Closing relationships can be determined using experimental data or CFD simulations. The equations of the porous body model [8] describe the averaged flow parameters that are continuously changing over the entire computational domain. The finite element method (FEM) [9] or the finite volume method (FVM) [10] applied for discretizing the computational domain allows for calculating geometries of arbitrary shape and does not require the use of structured grids, in contrast to the sub-channel approach. This approach allows for a determination of the average parameters at any point in the computational domain.

CFD models use a direct solution of NS equations (DNS approach) or their approximate solution using large eddy simulation (LES), Reynolds-averaged Navier–Stokes (RANS) approaches [11].

The DNS approach requires huge computational resources. It is expected that at the current rate of development of these resources, it will become available for solving problems of flow around a typical civil aircraft or car, with the number of grid nodes and the number of time steps at approximately 10¹⁶ and 10^7.7, respectively, no earlier than 2080 [12].

The LES approach, in which “large vortices” are described accurately, and approximate relations are used for “small” ones, requires noticeably less computational resources compared to the DNS approach [12]. However, the required resources remain very significant, so it is mainly used to validate the results obtained by the RANS approach, as well as to solve problems where the use of the RANS approach leads to incorrect results [12].

The RANS approach is based on the solution of the averaged NS equations and requires the use of semi-empirical turbulence models [12]. Its popularity is due to the relatively low cost of computing resources compared to using LES and even more so DNS approaches. Nevertheless, even using the RANS approach requires a huge amount of computational resources when simulating the core of a nuclear power plant [3]. RANS turbulence models are excellent for modeling fragments of an investigated object [13] to determine the required parameters or to determine the closing relationships of sub-channel models and a porous body model in the absence of experimental data. Also, at the moment, there is a problem with validating the results obtained using CFD codes on full-scale models, such as the reactor core [14].

In this paper, to describe thermo-hydraulic processes, an approach based on the model of an anisotropic porous body is used. The possibility of determining the closing relations of this model using a special integral turbulence model is shown. Calculated relationships are obtained for the parameters of the integral turbulence model for the main types of rod structures used in heat exchange equipment.

2. Anisotropic Porous Body Model

Hereinafter, the following notation is used:

$\mathsf{\Psi }$—variable; $\overline{\mathsf{\Psi }}$—time-averaged variable value; $\langle \mathsf{\Psi }\rangle$—volume-averaged variable value; ${\langle \mathsf{\Psi }\rangle }^i$—fluid volume-averaged variable value; and ${}^{\vee }\mathsf{\Psi }$—deviation from the value of the fluid volume-averaged variable. Averaging is carried out within the limits of a regular cell of the considered assembly of rods.

Korsun presented the equations of the anisotropic model of a porous body [8]. They consist of the equations of mass, momentum, and energy conservation and describe the motion of an averaged incompressible flow in a porous structure:

$$\frac{\partial }{\partial x_j}\varphi {\langle {\overline{u}}_j\rangle }^i=0$$

(1)

$$\frac{\partial }{\partial \tau }\rho \varphi {\langle \overline{u_i}\rangle }^i+\frac{\partial }{\partial x_j}\rho \varphi {\langle \overline{u_j}\overline{u_i}\rangle }^i=\rho \varphi g_i+{\overline{f}}_i+\frac{\partial }{\partial x_j}{\overline{T}}_{ij}^{\mathrm{eff}}$$

(2)

$$\rho c_p\varphi \frac{\partial {\langle \overline{t}\rangle }^i}{\partial \tau }+\frac{\partial }{\partial x_j}\rho c_p\varphi {\langle \overline{u_j}\overline{t}\rangle }^i=\frac{\partial }{\partial x_j}{\lambda }_{ij}^{\mathrm{eff}}\frac{\partial {\langle \overline{t}\rangle }^i}{\partial x_i}+\varphi \overline{q_V},$$

(3)

where: $u$—velocity; $\varphi$—porosity; $\rho$—density; $g_i$—components of the gravitational acceleration vector; ${\overline{f}}_i$—components of the vector of the bulk drag force; $T_{ij}^{\mathrm{eff}}$—components of the effective stress tensor; $c_p$—specific heat capacity; $t$—temperature; ${\lambda }_{ij}^{\mathrm{eff}}$—effective thermal conductivity tensor; and $q_V$—volumetric energy release in a liquid. The last term is also used to account for the heat gain from a solid to a liquid. Index $i=1,2,3$, and over repeated indices $j$ are summed from 1 to 3.

2.1. Closing Relations of the Anisotropic Model of a Porous Body

To determine the values of the averaged velocity, pressure, and temperature, the system of Equations (1)–(3) must be closed, namely, to define the form of the effective stress $T_{ij}^{\mathrm{eff}}$ and effective thermal conductivity ${\lambda }_{ij}^{\mathrm{eff}}$ tensors, and the volumetric resistance force $f_i$.

The bulk drag force is defined as follows [8]:

$${\overline{f}}_i=k_{ij}{\langle {\overline{u}}_j\rangle }^i,$$

(4)

where: $k_{ij}=k_{\xi \xi }{\delta }_{ij}+(k_{\eta \eta }-k_{\xi \xi })n_i{n}_j$—resistance tensor components; $k_{\eta \eta }(\beta )$, $k_{\xi \xi }(\beta )$—principal components of the tensor in the directions along and across the rods:

$$k_{\xi \xi }=\frac{c_{d\perp }}{2d_{\mathrm{h}}}\varphi \rho \left|{\langle \overline{\overrightarrow{u}}\rangle }^i\right|k_{\perp }(\beta )$$

(5)

$$k_{\eta \eta }=\frac{c_{d\uparrow }}{2d_{\mathrm{h}}}\varphi \rho \left|{\langle \overline{\overrightarrow{u}}\rangle }^i\right|k_{\uparrow }(\beta )$$

(6)

where: $c_{d\uparrow }(\mathrm{Re})$ and $c_{d\perp }(\mathrm{Re})$ are the coefficients of hydraulic resistance for longitudinal and transverse flow around the bundles of rods; $d_{\mathrm{h}}$—hydraulic diameter; $\left|\overrightarrow{u}\right|$—velocity module; and $k_{\perp }(\beta )$ and $k_{\uparrow }(\beta )$—correction factors.

Correction factors $k_{\perp }(\beta )$ and $k_{\uparrow }(\beta )$ were experimentally determined earlier and have the following form [15]:

$$k_{\perp }(\beta )={(1-\beta )}^{2.9},$$

(7)

$$k_{\uparrow }(\beta )={\beta }^{2.9}.$$

(8)

Here $\beta ={\cos }^2\phi$, where $\phi$ is the angle between the average velocity vector $\overrightarrow{u}$ and the unit vector $\overrightarrow{n}$.

For the components of the effective stress tensor, using the theory of matrix polynomials [16], Korsun obtained the following dependence [8]:

$$\begin{array}{l}{\overline{T}}_{ij}^{\mathrm{eff}}=-\varphi {\langle p{\delta }_{ij}\rangle }^i+\langle \mu (\frac{\partial {\langle {\overline{u}}_i\rangle }^i}{\partial x_j}+\frac{\partial {\langle {\overline{u}}_j\rangle }^i}{\partial x_i})\rangle -\varphi {\langle \rho \overline{{u^{\prime }}_i{u^{\prime }}_j}\rangle }^i-\varphi {\langle \rho \overline{{}^{\vee }{u}_i{}^{\vee }{u}_j}\rangle }^i=\\ =-\varphi {\langle P\rangle }^i+{\mu }_{\mathrm{eff}}(\frac{\partial {\langle {\overline{u}}_i\rangle }^i}{\partial x_j}+\frac{\partial {\langle {\overline{u}}_j\rangle }^i}{\partial x_i})\end{array}$$

(9)

$${\mu }_{\mathrm{eff}}=\left(\nu +{\langle {\nu }_t\rangle }^i\right)\left(\varphi -{\varphi }_1\right)\rho$$

(10)

$${\langle P\rangle }^i={\langle p{\delta }_{ij}\rangle }^i+\rho c_p{}^{\mathrm{t}}{\left|{\langle \overline{\overrightarrow{u}}\rangle }^i\right|}^2{\delta }_{ij},$$

(11)

where: $P$—effective pressure in the flow, consisting of thermodynamic pressure and pressure due to turbulent pulsations; ${\mu }_{\mathrm{eff}}$—effective dynamic viscosity; $-\varphi {\langle \rho \overline{{}^{\vee }{u}_i{}^{\vee }{u}_j}\rangle }^i$—is responsible for the transfer of impulse by the deflection velocities, which appears as a result of averaging the initial NS equations over the volume: $c_p{}^t$—pressure coefficient; ${\delta }_{ij}$—the unit tensor; $\nu$—kinematic viscosity; ${\nu }_t$—kinematic turbulent viscosity; and ${\varphi }_1$—porosity with dense packing of rods, when the spacing of the lattice of rods coincides with their diameter.

The components of the effective thermal conductivity tensor are determined by the following formula [17]:

$${\lambda }_{ij}^{\mathrm{eff}}=((\lambda +\rho c_p\frac{{\langle {\nu }_t\rangle }^i}{{\mathrm{Pr}}_t})(\varphi -{\varphi }_1)+{\lambda }^{\mathrm{rod}}){\delta }_{ij}+((\lambda +\rho c_p\frac{{\langle {\nu }_t\rangle }^i}{{\mathrm{Pr}}_t}){\varphi }_1-{\lambda }^{\mathrm{rod}})n_i{n}_j,$$

(12)

where: $\lambda$—thermal conductivity; ${\mathrm{Pr}}_t$—turbulent Prandtl number; ${\lambda }^{\mathrm{rod}}$—thermal conductivity through the rods; and $n_i$—the component of the unit vector coinciding with the axis of the bars.

The last two Equations for the components of the effective stress and effective thermal conductivity tensors are united by the need to determine the integral coefficient of effective kinematic viscosity ${\langle {\nu }_t\rangle }^i$. This coefficient can be determined using the integral turbulence model [18].

2.2. The State of Work on the Development of an Integral Turbulence Model

An analysis of the current state of work on the development of an integral turbulence model is presented in [19]. Most of the existing models have not been fully developed, and there is no information about testing and their practical application.

One of the most recent and sufficiently developed models is the model of Pedras and de Lemos [7]. This model uses the local k−ε turbulence model [20,21] as the base model. The authors integrate the local model over the volume in order to obtain an integral model of turbulence [7]. However, the authors of [7] use a number of unjustified assumptions when deriving the resulting equations for kinetic energy and dissipation rate of the kinetic energy of turbulence. Therefore, in the expression for the relationship between the averaged kinetic energy of turbulence and the rate of its dissipation with the averaged turbulent viscosity, the authors use a constant coefficient $c_{\mu }$, the value of which is taken from the usual k−ε turbulence model. To determine the closing relations of the developed model, the authors performed CFD calculations of only the transverse flow around the bundle of rods rectangular [22], ellipsoidal [23], and round [7] cross-section. The longitudinal flow around the rod bundles was not considered. The issues of anisotropy of fluid transfer and effective transfer coefficients in porous structures were also not considered. Examples of the practical application of the developed model were not presented. Due to the shortcomings of the integral model of Pedras and de Lemos [7], it was decided to develop a different integral turbulence model [18], also based on the k−ε turbulence model [10].

3. Integral k−ε Turbulence Model

The integral k−ε turbulence model [18] includes two equations for describing the transfer of the integral kinetic energy of turbulence k and its dissipation ε. The model was obtained by averaging over the space of the initial local equations of the standard k−ε turbulence model [10,24].

The final expressions for the equations of the integral k−ε turbulence model for rod assemblies are as follows:

$$\frac{\partial }{\partial \tau }\varphi {\langle k\rangle }^i+{\langle {\overline{u}}_j\rangle }^i\frac{\partial }{\partial x_j}\varphi {\langle k\rangle }^i=\frac{\partial }{\partial x_j}(-q_j^k)+P_u+P_d-\varphi {\langle \epsilon \rangle }^i$$

(13)

$$\frac{\partial }{\partial \tau }\varphi {\langle \epsilon \rangle }^i+{\langle {\overline{u}}_j\rangle }^i\frac{\partial }{\partial x_j}\varphi {\langle \epsilon \rangle }^i=\frac{\partial }{\partial x_j}(-q_j^{\epsilon })+C_{\epsilon }\frac{{\langle \epsilon \rangle }^i}{{\langle k\rangle }^i}\left[P_u+P_d-\varphi {\langle \epsilon \rangle }^i\right]$$

(14)

where: k—kinetic energy of turbulence; $q_j^k$—the flow of kinetic energy of turbulence; $P_u$—generation of turbulence due to averaged motion; $P_d$—generation of turbulence by the deflection velocities that appear as a result of averaging the initial local equations of the standard k−ε model of turbulence over the volume: $P_d=\langle -{}^{\vee }\overline{{u^{\prime }}_i{u^{\prime }}_j}\frac{\partial {}^{\vee }{\overline{u}}_i}{\partial x_j}\rangle$; $\epsilon$—the rate of dissipation of the kinetic energy of turbulence; $q_j^{\epsilon }$—flow rate of dissipation of kinetic energy of turbulence; and $C_{\epsilon }$—coefficient equal to 1.44.

The value of the integral kinematic turbulent viscosity can be determined by the following formula:

$${\langle {\nu }_t\rangle }^i={\langle c_{\nu }\rangle }^i\frac{{\langle k\rangle }^i{}^2}{{\langle \epsilon \rangle }^i}$$

(15)

where: ${\langle c_{\nu }\rangle }^i$ is the proportionality coefficient connecting the integral coefficient of the kinematic turbulent viscosity ${\langle {\nu }_t\rangle }^i$ with the average values ${\langle k\rangle }^i$ of the kinetic energy of turbulence and the rate of dissipation of the kinetic energy of turbulence ${\langle \epsilon \rangle }^i$. This coefficient is the closing relation in the model and needs to be determined. It should be noted that, in the standard k−ε turbulence model, this coefficient is constant and equal to 0.09 [10,24]. In the integral turbulence model, obtained by averaging the equations of the standard k−ε turbulence model, this coefficient does not have to be strictly constant and, therefore, its value must be refined.

The fluxes of turbulence kinetic energy ${\overrightarrow{q}}^k$ and the rate of its dissipation ${\overrightarrow{q}}^{\epsilon }$ are described by the following dependencies:

$${\overrightarrow{q}}^k=-(a_{mol}^k+a_t^k+a_{def}^k)\nabla {\langle k\rangle }^i=-a^k\nabla {\langle k\rangle }^i$$

(16)

$${\overrightarrow{q}}^{\epsilon }=-(a_{mol}^{\epsilon }+a_t^{\epsilon }+a_{def}^{\epsilon })\nabla {\langle \epsilon \rangle }^i=-a^{\epsilon }\nabla {\langle \epsilon \rangle }^i$$

(17)

where: $a_{mol}^k$—effective coefficient of transfer of kinetic energy due to molecular transfer; $a_t^k$—effective coefficient of transfer of kinetic energy due to turbulent transfer; $a_{def}^k$—effective coefficient of transfer of kinetic energy due to deflection velocities; $a_{mol}^{\epsilon }$—effective transfer coefficient of the rate of dissipation of kinetic energy due to molecular transfer; $a_t^{\epsilon }$—effective coefficient of transfer of the rate of dissipation of kinetic energy due to turbulent transfer; and $a_{def}^{\epsilon }$—effective coefficient of transfer of the rate of dissipation of kinetic energy due to the rates of deflection.

The effective transfer coefficients of kinetic energy $a^k$ and the rate of its dissipation $a^{\epsilon }$ are described by the deflection velocities:

$$a_{mol}^k=\nu \varphi$$

(18)

$$a_{{}_{mol,ij}}^{\epsilon }=\nu (\varphi -{\varphi }_1){\delta }_{ij}+\nu {\varphi }_1{n}_i{n}_j$$

(19)

$$a_t^k=\frac{{\langle {\nu }_t\rangle }^i}{{\mathrm{Pr}}_t^k}\varphi$$

(20)

$$a_{t,ij}^{\epsilon }=\frac{{\langle {\nu }_t\rangle }^i}{{\mathrm{Pr}}_t^{\epsilon }}\left[(\varphi -{\varphi }_1){\delta }_{ij}+{\varphi }_1{n}_i{n}_j\right]$$

(21)

where: $n_i$—component of the unit vector, coinciding with the axis of the rods; and ${\mathrm{Pr}}_t^k$, ${\mathrm{Pr}}_t^{\epsilon }$—the turbulent Prandtl numbers.

Generations of turbulence by averaged motion $P_u$ and by deflection velocities $P_d$ are described by the following dependencies:

$$P_u=\varphi {\langle -\overline{{u^{\prime }}_i{u^{\prime }}_j}\rangle }^i\frac{\partial {\langle {\overline{u}}_i\rangle }^i}{\partial x_j}$$

(22)

$$P_d=\langle -{}^{\vee }\overline{{u^{\prime }}_i{u^{\prime }}_j}\frac{\partial {}^{\vee }{\overline{u}}_i}{\partial x_j}\rangle$$

(23)

$$P_d=P_{d\uparrow }·\beta +P_{d\perp }(1-\beta )$$

(24)

$$P_{d\uparrow }=\frac{N}{\rho }·{\langle {\alpha }_v\rangle }^i$$

(25)

$$P_{d\perp }=0.28{\varphi }^2\frac{{\langle k\rangle }^i·\left|{\langle \overline{\overrightarrow{u}}\rangle }^i\right|}{\sqrt{K}}$$

(26)

where: $P_{d\uparrow }$—the generation of turbulence by the deflection velocities in the longitudinal flow; $P_{d\perp }$—generation of turbulence by deflection velocities in the transverse flow [7]; and ${\langle {\alpha }_v\rangle }^i$—the share of energy transferred to turbulent pulsations. This factor can be determined in the course of computational modeling. The variable $N$—the total power equal to the ratio of the work of pushing the coolant through the assembly per unit liquid volume; $K$—permeability of medium.

Using Equation (25) and considering that the share of power transmitted to turbulent pulsations $N_t$ is determined by the rate of dissipation of kinetic turbulent energy ε, the following formula was obtained to determine the value of the integral coefficient ${\langle {\alpha }_v\rangle }^i$:

$${\langle {\alpha }_v\rangle }^i=\frac{N_t}{N}=\frac{\rho {\langle \epsilon \rangle }^{i}V}{N}$$

(27)

where: $V$—the volume of the regular cell.

4. Determination of the Closing Relations for the Integral Turbulence Model

To close the system of equations for the anisotropic model of a porous body (1)–(3), it is necessary to determine the values of the components of the effective stress tensor ${\overline{T}}_{ij}^{\mathrm{eff}}$ and the effective thermal conductivity tensor ${\lambda }_{ij}^{\mathrm{eff}}$. These coefficients can be calculated through the integral coefficient of kinematic turbulent viscosity ${\langle {\nu }_t\rangle }^i$, which is determined using the integral turbulence model (13)–(15), for the closure of which it is necessary to determine the coefficients ${\langle c_{\nu }\rangle }^i$ and ${\langle {\alpha }_v\rangle }^i$.

One can determine the value of these coefficients experimentally or with help of CFD calculation [17]. In this work, the coefficients ${\langle c_{\nu }\rangle }^i$ and ${\langle {\alpha }_v\rangle }^i$ were determined using systematic computational studies in a wide range of Reynolds numbers and porosities of a rod assembly. The cases of longitudinal and transverse flow around bundles of rods were considered. The calculations were carried out using the standard local k−ε turbulence model [24] in the ANSYS CFX [10]. ANSYS CFX is a general-purpose CFD package for simulating fluid and gas flow considering turbulence, heat transfer, interfacial interactions, chemical reactions, and combustion.

Closing relations for the integral turbulence model were obtained for smooth bundles of rods, fuel assemblies with a spacer grid, or wire-wrapped fuel bundles.

4.1. Computational Studies of the Flow around Smooth Bundles of Rods

The transverse and longitudinal flow of water with a temperature of approximately 25 °C was simulated for the assembly of round rods (Figure 1 and Figure 2). The diameter of the rods was 10 mm. Reynolds number varied from 10⁴ to 10⁶. The ordinary definition of the Reynolds number in a bundle of rods was used. For a longitudinal flow, the Reynolds number was calculated by the average velocity in the cross-section of the rod bundle. In this case, the hydraulic diameter of an elementary regular bundle cell was used as a characteristic size. The Reynolds number was calculated by the average velocity in the narrow bundle section and the rod diameter in the case of a transverse flow. The porosity value varied from 0.4 to 0.8 by changing the pitch value. The length of the simulated section in the case of longitudinal flow was 20 mm.

The calculations were carried out using the standard local k−ε turbulence model [24] in the ANSYS CFX code [10]. The scalable wall functions were used in the calculation [10]. The scheme of the second order of accuracy for the convective term in the equations of the law of conservation of momentum and the second order of accuracy for the convective term in the equation of turbulent characteristics were used. Unbalance values less than 1% and maximum local residuals not more than 10⁻⁵ were chosen as criteria of convergence.

Periodic boundary conditions were set at the end boundaries of the computational domain with the specified flow rate. The boundary symmetry condition was set on the lateral surfaces. A no-slip condition was set on the wall.

The computational mesh was created using the ANSYS MESH program [25,26]. The size of the grid element was set in such a way as to obtain the grid convergence of the solution for the required integral parameters ${\langle {\nu }_t\rangle }^i$, ${\langle k\rangle }^i$, ${\langle \epsilon \rangle }^i$. Additionally, the pressure drop on the simulated domain was controlled. This drop was used to check simulation results by comparing it with the calculation based on empirical dependences from the reference book [27]. The resulting deviations of the pressure drop from those calculated from the experimental dependences did not exceed 10% for longitudinal flow and 20% for transverse flow, which are within the error of empirical dependences.

For the case of longitudinal flow, the number of mesh elements ranged from 353,880 to 788,000 for a triangular bundle of rods and from 528,560 to 1,163,640 for a square bundle of rods, depending on the Re number and porosity. In the direction of the fluid flow, the region was divided into 40 elements, with decreasing the mesh element size towards the boundaries (inlet and outlet). The y+ value varied in the range from 1 to 180.

For the case of transverse flow, the number of mesh elements ranged from 240,856 to 492,689 for a triangular bundle of rods and from 101,357 to 630,687 for a square bundle of rods, depending on the Re number and porosity. The y+ value varied in the range from 0.3 to 150.

The stationary distribution of velocities, pressure, kinetic energy of turbulence, and the rate of dissipation of kinetic energy of turbulence within the calculated volume were obtained. The integral characteristics of the flow were then obtained for cell-averaged values of the kinetic energy of turbulence k, the rate of dissipation of kinetic energy of turbulence ε, the coefficient of kinematic turbulent viscosity ${\nu }_t$, and integral coefficient ${\langle c_{\nu }\rangle }^i$ in expression (15).

The dependences of the integral coefficient ${\langle c_{\nu }\rangle }^i$ on the Re number for longitudinal and transverse flow around the assembly of round rods are shown in Figure 3 and Figure 4.

The integral coefficient ${\langle c_{\nu }\rangle }^i$ is a function of the Reynolds number for longitudinal flow around the assembly of round rods:

$${\langle c_{\nu }\rangle }^i=0.0134·{\left({\langle \mathrm{Re}\rangle }^i\right)}^{0.25}.$$

(28)

In Equation (28), Re is the Reynolds number.

The integral coefficient ${\langle c_{\nu }\rangle }^i$ is practically constant for a transverse flow around the assembly of round rods and can be taken equal to 0.117.

The dependence of the integral coefficient ${\langle c_{\nu }\rangle }^i$ on the Re number for an arbitrary angle of flow around the assembly, in the first approximation, can be represented as follows:

$${\langle c_{\nu }\rangle }^i=0.0134·{\left({\langle \mathrm{Re}\rangle }^i\right)}^{0.25}·\beta +0.117·(1-\beta ).$$

(29)

In the local k−ε turbulence model, the coefficient $c_{\nu }$ is constant and equals 0.09 [28]. In the integral turbulence model, this coefficient ${\langle c_{\nu }\rangle }^i$ is not a constant and depends both on the flow direction and on the Reynolds number but does not depend on the porosity of the structure or the type of rod assembly.

Figure 5 shows the dependence of the integral coefficient ${\langle {\alpha }_v\rangle }^i$, on the Re number and the porosity of the assembly, which was calculated using Equation (27).

The figure shows that the integral coefficient ${\langle {\alpha }_v\rangle }^i$ increases from 0.2 to 0.7 with an increase in the Re number when it changes in the range from 10⁴ to 10⁶. It should be noted that the porosity of the structure and the type of packing of the rods have little effect on the value of ${\langle {\alpha }_v\rangle }^i$. The results can be approximated by the following relationship:

$${\langle {\alpha }_v\rangle }^i=0.1·\ln \left({\langle \mathrm{Re}\rangle }^i\right)-0.63.$$

(30)

The set of closing relations for smooth bundles of rods is determined by Equation (29) for the integral coefficient ${\langle c_{\nu }\rangle }^i$ and Equation (30) for the integral coefficient ${\langle {\alpha }_v\rangle }^i$.

4.2. Computational Studies of the Flow around Fuel Assemblies with a Spacer Grid and Wire-Wrapped Fuel Bundles

In fuel assemblies, fuel rod spacing is carried out using special spacer elements, which are spacer grids or wire winding on fuel rods. The spacer elements affect the flow of the coolant in the fuel assembly and require a corresponding correction of the integral turbulence model.

4.2.1. The Set of Closing Relations for Fuel Assemblies with a Spacer Grid

The presence of the spacer grid leads to a change in the form of the resistance tensor and the generation of turbulence in the spacer grid. The set of closing relations outside of the spacer grid area is the same as the set of closing relations for smooth bundles of rods (Equations (29) and (30)).

In the area of the spacer grid, the values of the tensor of resistance and the generation of turbulence change. The coefficient of hydraulic resistance for the spacer grid can be calculated using the following formula [29]:

$${\xi }_D=\frac{(1-m)(1.4-0.5m)}{m^2}+0.11{\left(\frac{64}{{\mathrm{Re}}_p}+\frac{\Delta }{d_{\mathrm{h}}}\right)}^{0.25}\frac{l_p}{m^2{d}_{\mathrm{h}}} ,$$

(31)

In Equation (31), ${\xi }_D$ is the constraint coefficient and equals the ratio of the cross-sectional area to the area of the flow area of the channel; $m=1-{\xi }_D$ is the narrowing coefficient; l_p is the height of the spacer grid; and Δ is the roughness of the spacer grid.

The rate of turbulent production due to deflection velocities $P_d$ can be calculated using the following formula:

$$P_D(z)={\langle {\alpha }_v\rangle }^i{\xi }_{Di}\frac{\rho {\left({\langle {\overline{u}}_i\rangle }^i\right)}^2}{2}\frac{{\langle {\overline{u}}_i\rangle }^i}{l_p}\varphi \equiv P_{Di}\frac{1}{l_p} , z_i-\frac{l_p}{2}\le z\le z_i+\frac{l_p}{2},$$

(32)

In Equation (32), $P_{Di}\equiv {\langle {\alpha }_v\rangle }^i{\xi }_{Di}\frac{\rho {\left({\langle {\overline{u}}_i\rangle }^i\right)}^2}{2}{\langle {\overline{u}}_i\rangle }^i\varphi$ is the total turbulence energy generated by the spacer grid; $P_D(z)$ is not equal to zero only in the spacer grid region; ${\xi }_{Di}$ is the longitudinal coefficient of hydraulic resistance for a spacer grid; ${\langle {\overline{u}}_i\rangle }^i$ is the inlet flow velocity in the spacer grid region; and l_p is the height of the spacer grid. The value of the integral coefficient ${\langle {\alpha }_v\rangle }^i$ must be determined for a spacer grid with the help of CFD modeling.

4.2.2. The Set of Closing Relations for Wire-Wrapped Fuel Bundles

In the area of the wire-wrapped fuel bundles, the values of the tensor of resistance change. The coefficients of hydraulic resistance for the wire-wrapped fuel bundles can be calculated using the following formulas [30,31] for longitudinal and transverse flow:

$$c_{\mathrm{d}\uparrow }=c_0\left(1+600{\left(\frac{d}{T}\right)}^2\left(\frac{s}{d}-1\right)\right)$$

(33)

$$c_{\mathrm{d}\perp }=c_0\left[\left(1+600{\left(\frac{d}{T}\right)}^2\left(\frac{s}{d}-1\right)\right)\left(\frac{0.16}{{(\raisebox{1ex}{$s$}\!\left/ \!\raisebox{-1ex}{$d$}\right.-1)}^2}+24·\frac{s}{d}\right)\right],$$

(34)

In Equations (33) and (34), $c_0$ is the coefficient of longitudinal hydraulic resistance for smooth rods; T is the wire lead length; d is the rod diameter; and s is the rod pitch.

It was shown earlier [32] that to describe the integral coefficients ${\langle c_{\nu }\rangle }^i$ and ${\langle {\alpha }_v\rangle }^i$ one can use the dependence for smooth bundles of rods for wire-wrapped fuel bundles.

5. Numerical Implementation and Testing of the Developed Model

For the numerical solution of the equations of the porous body model and the integral turbulence model, FEM is used. The FEM is implemented using the DOLFIN library of the FEniCS package [33].

5.1. Variational Formulation of the Problem

The so-called weak form of the boundary value problem for the equations of the porous body model (1)–(3) together with the integral turbulence model (13)–(15) was obtained by multiplying the original differential equations by test functions and then using integration by parts [9] (to simplify the notation, the averaging signs were omitted):

$$\begin{array}{l}F((u,p),({\upsilon }_u,{\upsilon }_p))=(1/d\tau ){\displaystyle \underset{\mathsf{\Omega }}{\int }(u-u_0)}·{\upsilon }_u\mathrm{d}x+{\displaystyle \underset{\mathsf{\Omega }}{\int }(u·\nabla u)}·{\upsilon }_u\mathrm{d}x\\ +{\nu }_{eff}{\displaystyle \underset{\mathsf{\Omega }}{\int }(\nabla u+\nabla u^{\mathrm{T}})}·\nabla {\upsilon }_u\mathrm{d}x-\frac{1}{\rho }{\displaystyle \underset{\mathsf{\Omega }}{\int }p}·div{\upsilon }_u\mathrm{d}x\\ +{\displaystyle \underset{\mathsf{\Omega }}{\int }divu·{\upsilon }_p}\mathrm{d}x-\frac{1}{\rho \varphi }{\displaystyle \underset{\mathsf{\Omega }}{\int }f}·{\upsilon }_u\mathrm{d}x-{\displaystyle \underset{\mathsf{\Omega }}{\int }g}·{\upsilon }_u\mathrm{d}x\\ -{\nu }_{eff}{\displaystyle \underset{\partial {\mathsf{\Omega }}_N}{\int }(\nabla u^{\mathrm{T}}·n)}·{\upsilon }_u{\mathrm{ds}}_{outlet}+\frac{1}{\rho }{\displaystyle \underset{\partial {\mathsf{\Omega }}_N}{\int }(p_{outlet}·n)·{\upsilon }_u}{\mathrm{ds}}_{outlet}=0 \end{array}$$

(35)

$$\begin{array}{l}F(k,{\upsilon }_k)=(1/d\tau ){\displaystyle \underset{\mathsf{\Omega }}{\int }(k-k_0)}·{\upsilon }_k\mathrm{d}x+{\displaystyle \underset{\mathsf{\Omega }}{\int }(u·\nabla k)}·{\upsilon }_k\mathrm{d}x+a_k{\displaystyle \underset{\mathsf{\Omega }}{\int }\nabla k}·\nabla {\upsilon }_k\mathrm{d}x\\ -{\displaystyle \underset{\mathsf{\Omega }}{\int }\frac{1}{\varphi }·P_u·{\upsilon }_k}\mathrm{d}x-{\displaystyle \underset{\mathsf{\Omega }}{\int }\frac{1}{\varphi }·P_d}·{\upsilon }_k\mathrm{d}x+{\displaystyle \underset{\mathsf{\Omega }}{\int }\epsilon }·{\upsilon }_k\mathrm{d}x=0 \end{array}$$

(36)

$$\begin{array}{l}F(\epsilon ,{\upsilon }_{\epsilon })=(1/d\tau ){\displaystyle \underset{\mathsf{\Omega }}{\int }(\epsilon -{\epsilon }_0)}·{\upsilon }_{\epsilon }\mathrm{d}x+{\displaystyle \underset{\mathsf{\Omega }}{\int }(u·\nabla \epsilon )}·{\upsilon }_{\epsilon }\mathrm{d}x+{\displaystyle \underset{\mathsf{\Omega }}{\int }(a_{\epsilon }}·\nabla \epsilon )·\nabla {\upsilon }_{\epsilon }\mathrm{d}x\\ -{\displaystyle \underset{\mathsf{\Omega }}{\int }{C}_{\epsilon }\frac{1}{k}\left[\frac{1}{\varphi }·P_u·\epsilon +\frac{1}{\varphi }·P_d·\epsilon -\epsilon ·\epsilon \right]}·{\upsilon }_{\epsilon }\mathrm{d}x=0 \end{array}$$

(37)

$$\begin{array}{l}F(t,{\upsilon }_t)=(1/d\tau )\rho c_p{\displaystyle \underset{\mathsf{\Omega }}{\int }(t-t_0)}·{\upsilon }_t\mathrm{d}x+\rho c_p{\displaystyle \underset{\mathsf{\Omega }}{\int }(u·\nabla t)}·{\upsilon }_t\mathrm{d}x\\ +{\displaystyle \underset{\mathsf{\Omega }}{\int }({\lambda }_{eff}}·\nabla t)·\nabla {\upsilon }_t\mathrm{d}x-{\displaystyle \underset{\mathsf{\Omega }}{\int }{q}_V}·{\upsilon }_t\mathrm{d}x=0,\end{array}$$

(38)

In Equations (35)–(38), ${\upsilon }_u, {\upsilon }_p, {\upsilon }_k, {\upsilon }_{\epsilon }, {\upsilon }_t$ are test functions; $\mathsf{\Omega }$ is the area; $\partial {\mathsf{\Omega }}_N$ is the area boundary; n is the normal to the surface; ${\nu }_{eff}$ is the turbulent viscosity.

Equation (35) is written for the case of using a “coupled” approach to solving the equation of the law of conservation of momentum, which solves all equations for all desired variables for a particular element and repeats this process for all elements into which the computational domain is divided.

Taylor–Hood finite elements P2-P1 were used to discretize the equation of the law of conservation of momentum [34]. Lagrange finite elements were used for the equations of transfer of the kinetic energy of turbulence (36), the rate of dissipation of the kinetic energy of turbulence (37), and energy (38).

5.2. Scheme for Solving a System of Equations

At each time step, a nonstationary problem is solved with a given number of iterations. At each internal iteration, the equations for the velocity and pressure are first solved. The equations for the turbulence parameters are then solved. Next, the equations for the conservation of energy are solved. The solution of the equations occurs until the specified maximum number of iterations for the corresponding equation is exceeded, or upon reaching the specified calculation accuracy, depending on which comes first. Additionally, the imbalance value is calculated for the equation of the law of conservation of momentum. At the end of the iterative process, the results are saved in data file format provides by open-source, freely available Visualization Toolkit (VTK) [35].

5.3. Numerical Model Testing on a 19-Bar Experimental Assembly

Simulation of the experiment was performed. The experiment was previously carried out at the Institute for Physics and Power Engineering (IPPE) on an experimental 19-bar assembly [36]. In the experiment, Na−K was used as a coolant. The hexagonal fuel assembly used in the experiment consisted of 19 unheated smooth rods embedded in a hexagonal channel. The rod diameter was 19 mm and its relative pitch was 1.17. The apothem of the hexagonal channel was 49.4 mm. The rods were spaced by the upper and lower gratings [36]. The distance between gratings was 745 mm. The lower ends of the rods were fixed in the spacer grid. The upper ends of the rods were extended beyond the working length of the model and sealed in the flange [36]. Blocking of the flow area was carried out using a plate. The plate was installed in the central part of the fuel assembly cross-section at a distance of 380 mm from the lower grate and overlapped the cross-section of the assembly to the middle of the peripheral rows of rods, forming a 55% central blockade of the flow area. The blockade thickness was 8 mm. The longitudinal velocity was measured along the length by the electromagnetic method. In the experiment, the values of the emf of the electromagnetic sensor were obtained for the central, intermediate, lateral, and angular channels (see Figure 6).

Stationary isothermal coolant flow around the assembly was simulated. The physical properties of the Na-K coolant were calculated at a temperature of 100 °C. The Re number was 3.3 × 10⁴. The calculation was carried out until the stationary distribution of all the required parameters was obtained. The time step value varied in the range 10⁻³–10⁻⁴ depending on the magnitude of the current imbalances for the equation of the law of conservation of momentum.

At the input, a boundary condition of the first type was set for the velocity, energy, kinetic energy of turbulence, and its dissipation rate. The values for the kinetic energy of turbulence and its dissipation rate were calculated through the value of the turbulence intensity, which was set to 5%. At the output, a boundary condition of the first type was set for the pressure and a zero gradient for all other parameters. The outlet pressure was set equal to 1 atm. The slip boundary condition for the velocity and zero gradients for all other parameters were set on the wall.

The construction of the geometry and mesh of the computational domain was carried out in the GMSH program [37]. The geometry of the computational domain is a hexagonal prism with a height of 745 and a wrench size of 98.8 mm. Mesh elements are tetrahedrons. The mesh consists of 6795 nodes and 37,224 elements. The mesh was divided into several sub-areas to set the sizes of mesh elements more accurately in the blockade area and behind it. The blocking area was simulated by high hydraulic resistance.

As a result of the calculation, the stationary distribution of velocities, pressure, kinematic turbulent viscosity, kinetic energy of turbulence, and dissipation rate of kinetic energy of turbulence were measured. The analysis of the results was carried out using the ParaView program [38].

Figure 7 shows a comparison of the results of calculating the longitudinal velocity along with the height with experimental data for the central element. Because for the central element there is no unevenness of the velocity along its perimeter, the authors of [36] give as experimental data the values of the longitudinal velocity averaged over all points 1–7 in Figure 8. The calculated values on the axis passing through the geometric center of the model are used as the simulation results.

When comparing the results of calculating the longitudinal velocity along the height with the experimental data for the intermediate element, the presence of uneven velocity along the perimeter was considered. The values at the points shown in Figure 8 were taken as the calculation results.

Figure 9 shows a comparison of the results of calculations of the longitudinal velocity along with the height with the experimental data at point 3 for the element.

The simulation results qualitatively agree with the experimental results. Behind the blockade, there is a vortex region with a reverse flow (Figure 7 and Figure 9). The maximum reverse flow velocity is on the bundle axis. The vortex region size is determined by the distance from the blockade to the rear stagnation point with zero velocity. Downstream of this region, the flow velocity gradually returns to the inlet velocity.

The calculated velocity values in the section before the blockade are in good agreement with the experimental data. The error in determining the maximum reverse flow velocity lies within 10–15%. At the same time, the calculated size of the vortex region differs from the experimental one by more than 40%. In addition, the calculated flow velocity returns to the initial value faster.

One of the reasons for the decrease in the vortex region size can be the blocking model used. The experimental blockade was made by a special plate, so it was completely impenetrable. In the porous body model used in the calculation, the blockade is simulated by setting a high hydraulic resistance in the blocking area, and thus it is partially permeable. The blockade permeability is approximately 3–5%. The presence of a leak flow reduces the size of the recirculation area. At 55% central blockage, a leakage flow of approximately 30% eliminates the reverse flow region [39].

The faster return of the calculated flow velocity to the initial state is possibly related to the boundary conditions at the outlet of the computational domain. In the experiment, a grating was located at the outlet of the bundle. It undoubtedly interfered with the transverse velocity equalization. The grating was not considered in the calculation and the boundary condition of constant pressure was used.

A detailed analysis of the reasons for the obtained discrepancies will be performed in further studies.

6. Conclusions

A method for calculating the closing relations of a porous body model using an integral turbulence model has been developed. During computational studies, the parameters of the integral turbulence model and the closing relations of the porous body model for assemblies of smooth rods were determined. The integral coefficient ${\langle c_{\nu }\rangle }^i$ is practically constant for a transverse flow around the assembly of round rods and can be taken equal to 0.117. The integral coefficient ${\langle c_{\nu }\rangle }^i$ is a function of the Reynolds number for longitudinal flow around the assembly of round rods. The integral coefficient ${\langle {\alpha }_v\rangle }^i$ increases from 0.2 to 0.7 with increasing Reynolds from 10⁴ to 10⁶ and weakly depends on the porosity of the structure and the type of packing of the rods. The presence of the spacer grid leads to a change in the form of the resistance tensor and the generation of turbulence in the spacer grid. The value of the integral coefficient ${\langle {\alpha }_v\rangle }^i$ must be determined for a spacer grid with the help of CFD modeling. The presence of wire-wrapped fuel bundles leads to a change in the form of the resistance tensor. The integral coefficients ${\langle c_{\nu }\rangle }^i$ and ${\langle {\alpha }_v\rangle }^i$ or wire-wrapped fuel bundles can be determined from the dependencies for smooth bundles of rods. The experiment simulation performed using the porous body model with the integral turbulence model confirmed the possibility of using this approach to describe thermal-hydraulic processes during fluid flow in rod structures. The scatter of the calculated values of coefficients ${\langle c_{\nu }\rangle }^i$ and ${\langle {\alpha }_v\rangle }^i$ may be attributed to the dependence of this coefficient not only on the Reynolds number, but also on the initial fluctuations in the flow and near the wall. Therefore, improving the model is possible, considering the stochastic relations according to [40,41,42,43,44,45,46].