Investigation of Point-Contact Strategies for CFD Simulations of Pebble-Bed Reactor Cores

Abstract

This study numerically investigated the effects of various contact strategies on the thermal hydraulic behavior within a structured bed of 100 explicitly modeled pebbles. Four contact strategies and two thermal hydraulic conditions were considered. The strategies to avoid contact singularities include decreasing the pebble diameter, increasing the pebble diameter, bridging the pebble surfaces near the contact region, and capping the pebble surfaces near the contact region. One strategy, Strategy 3a, which involves bridging with a cylinder equal to 10% of the pebble diameter, was selected as the baseline strategy because it addressed the contact singularity while minimizing the geometric changes that affect the bed porosity. The two thermal hydraulic conditions were full-power operation (Case 1) and pressurized loss of forced cooling or PLOFC (Case 2). Simulations of the conjugate heat transfer within the structured bed were performed using the Reynolds-averaged Navier–Stokes approach with the realizable k-$ϵ$ turbulence model and two-layer all $y^{+}$ wall treatment. The thermal-fluid quantities of interest were compared between the contact strategies for each case. In Case 1, the hydraulic behavior was sensitive to the contact strategy, with large differences in the pressure drop (30%) and volume-average velocity (4%). The thermal behavior was not sensitive, with less than a 0.5% difference across the strategies. To better understand the separate effects of each heat transfer mode, Case 2 was divided into the following subcases: conduction (Case 2a); conduction/radiation (Case 2b); and conduction/radiation/convection (Case 2c). Case 2a represents an early phase of the PLOFC transient. Case 2b represents an intermediate phase of the PLOFC transient, with the pebble temperatures sufficiently high for the radiative heat transfer to be non-negligible. Case 2c represents a late phase of the PLOFC transient after the establishment of the natural circulation of the heat transfer fluid. For Case 2, large differences in the contact strategy were observed only in Case 2a with only conduction. The difference in the maximum pebble temperature was 23% in Case 2a, 2% in Case 2b, and 0.3% in Case 2c.

1. Introduction

The primary objective of this study was to identify the optimal strategy to address the point-contact singularity commonly encountered when explicitly modeling pebble-bed reactor cores for computational fluid dynamics (CFD) simulations. Pebbles modeled as perfect spheres within structured or unstructured beds are in point-contact configurations with the adjacent pebbles (particle–particle) as well as internal reactor vessel and neutron reflector surfaces (particle–wall). This effort investigated structured beds.

While unstructured bed flow fields are different, the contact strategy options are the same regardless of if the beds are structured or unstructured. The most desirable contact strategies are those that prioritize minimizing the porosity error as much as possible while providing a thermal connection to enable specifying an effective point-contact thermal resistance. Therefore, the recommended contact strategies and sizing selection are expected to hold for unstructured beds. Also, a study investigating externally heated unstructured beds reached similar conclusions as this effort with internally heated pebbles regarding the opportunity cost of each strategy and the sizing recommendations [1].

However, the overall conclusions are expected to hold for unstructured beds. Such point singularities must be avoided to discretize the computational domain into finite-volume elements. The commonly used techniques to avoid these point contacts can be placed into two groups, which either decrease or increase the pebble volume. The methods that decrease the volume include decreasing the pebble diameter or chamfering the pebble near the point contact. Decreasing the pebble diameter is the most straightforward approach for model building, the least difficult to mesh, and often the preferred starting point. The first 3D explicitly modeled pebble bed with simultaneous heat transfer and fluid flow was investigated in 2004 via a pebble diametral reduction [2]. At that time, and today, the experimental techniques lack the ability to provide detailed local heat transfer measurements for contact strategy validation. The quasi-direct numerical simulation (DNS) results of a single face-centered cubic (FCC) unit cell were generated for model development and validation [3]. These quasi-DNS results used the diametral reduction strategy to avoid the point-contact singularity. The authors mentioned that their objective was to keep the inter-pebble gap reasonably small while adhering to the DNS mesh constraints. A 4.2% diametral reduction was chosen, a significantly greater decrease than the values studied in this paper (0.1–1%). The diametral reduction strategy is still used today in modern large-eddy simulation (LES) approaches with thousands of explicitly modeled pebbles [4] and DNS simulations with up to 150 pebbles [5]. Such a global geometric change was recognized to significantly alter the bed porosity and pressure drop [1]. Therefore, local geometric modifications are also of interest. The capping method involves chamfering the contacting pebbles and filling the chamfered volume with fluid [6].

There are also global and local increase volume methods that include increasing the pebble diameter or chamfering the contacting pebbles near the point contact and filling the nearby volume with a bridge of solid pebble material. Eleven explicitly modeled pebble beds were studied using global diametral increases [7,8]. However, a range of values did not consider estimating the sensitivity of the diametral increase and making suggestions about the desired amount of overlap. Two disadvantages of the diametral increase strategy include (1) having the potential to introduce new contacts with nearby adjacent pebbles, and (2) the fluid domain meshing is still complex in the small interstices near the original point-contact singularity [1].

Each method has a unique impact on the thermal hydraulic behavior of the bed. Therefore, the optimal contact strategy may not be the same for all the reactor conditions where the dominant mode of heat transfer changes from forced convection (steady-state full-power operation) to either conduction- or radiation-dominant. In this study, point-contact strategies were investigated for full-power operation at a typical pebble power density and for a pressurized loss of forced cooling (PLOFC) with a reactor trip event. Each contact strategy is illustrated in Figure 1.

• Strategy 1 (a-b-c): Decreasing the pebble diameter by $0.1\%$, $0.5\%$, and $1\%$. This is a simple model-building approach of slightly decreasing the pebble diameter. Decreasing the pebble by $0.1\%$ was the most computationally expensive approach due to the very small gaps between the pebbles.
• Strategy 2 (a-b-c): Increasing the pebble diameter by $0.1\%$, $0.5\%$, and $1\%$.
• Strategy 3 (a-b-c): Bridging—that is, chamfering contact pebbles with a given diameter ($10\%$, $15\%$, and $20\%$ in this study) and filling the volume with pebble material. This approach was considered the most similar to the true bed configuration when comparing its porosity to the true value, and thus it was selected as the baseline for comparison.
• Strategy 4 (a-b-c): Capping—that is, chamfering contact pebbles with a given diameter ($10\%$, $15\%$, and $20\%$ in this study).

Two reactor operating conditions were investigated for each contact strategy: (Case 1) full power and (Case 2) PLOFC. The Case 1 power density of $1\times {10}^7$ $\mathrm{W}/{\mathrm{m}}^{-3}$ is prototypical of full-power operation for tristructural isotropic (TRISO) fuel forms using high-assay low-enriched uranium and helium gas as the heat transport fluid [9]. The associated full-power values were specified as boundary conditions, including a 3 m/s⁻¹ fluid inlet velocity, a 250 °C inlet temperature, and an outlet pressure of 6 $\mathrm{M}\mathrm{Pa}$. Case 2 roughly corresponds to the PLOFC with a reactor trip event, and this case was motivated by an objective to better understand the contact strategy under the separate effects of each heat transfer mode of conduction, convection, and radiation. For this case, the pebbles were assumed to have a temporally constant power density equal to $7\%$ of the full-power density, which corresponds to approximately 20 s after the reactor trip [10].

Instead of temporally resolving the PLOFC transient, the event was divided into three unique and steady subcases that correspond to the distinct scenarios of the event. Case 2a represented a phase of the PLOFC transient when the heat transport fluid reaches a near-zero velocity during the expected flow reversal phenomena after circulator pump coastdown and control rod scram [11]. Case 2a only considered conduction as a means of energy transport to study how different the temperature and heat flux fields would be between contact strategies that result in pebble overlaps (Strategies 2 and 3) or pebble gaps (Strategies 1 and 4). In Case 2a, the pebble surface temperatures reached sufficiently high temperatures that would make thermal radiation a non-negligible mode of energy transport. Therefore, Case 2b included radiative heat transfer. Case 2c represented a late phase of the PLOFC transient after the establishment of the buoyancy-driven natural circulation of the heat transfer fluid. This case captured the contact strategy behavior when conduction, convection, and radiation occur simultaneously.

In this work, the pore-scale hydraulic diameter ($d_h$) and pore Reynolds number ($R{e}_p$) are defined as

$$d_h=d_p{\displaystyle \frac{{ϵ}_i}{1-{ϵ}_i}},$$

(1)

$$R{e}_p={\displaystyle \frac{{\rho }_f{U}_{int}{d}_h}{{\mu }_f}},$$

(2)

where ${ϵ}_i$ is the calculated porosity for each contact strategy, $d_p$ is the pebble diameter, and ${\rho }_f$ and ${\mu }_f$ are the density and dynamic viscosity of the heat transfer fluid (helium gas), respectively. $U_{int}$ is the average interstitial velocity between the pores related to the superficial velocity, U, through $U_{int}=U/ϵ$. The superficial velocity or the Darcian flow is defined as $U={\displaystyle \frac{Q}{A}}$, where Q is the volumetric flow rate, and A is the cross-sectional area of an empty bed.

The heat transfer coefficient $h_{sf}$ and the Nusselt number $N{u}_{sf}$ are defined as

$$h_{sf}={\displaystyle \frac{q_p^{″}}{\left({\overline{T}}_p-{\overline{T}}_f\right)}},$$

(3)

$$N{u}_{sf}={\displaystyle \frac{h_{sf}{d}_p}{{\lambda }_f}},$$

(4)

where $q_p^{″}$ is the pebble surface heat flux, and ${\overline{T}}_p$ and ${\overline{T}}_f$ are the volume-averaged pebble surface and fluid temperature. The fluid thermal conductivity is denoted as ${\lambda }_f$.

For each contact strategy, the calculated porosity, ${ϵ}_i$, was compared to an analytically computed true porosity using perfectly spherical pebbles with point singularity contact.

Table 1. Calculated porosity, pore-scale hydraulic diameter $d_h$, and porosity error relative to the true porosity for four different pebble contact strategies of various sizes.

Strategies	Pebble Volume [$\times {10}^{-3}$ m3]	Calculated Porosity	${\mathit{d}}_{\mathit{h}}$ [m]	Porosity Error [%]
1a—Decrease $0.1\%$	4.5678	0.4780	0.0458	0.3284
1b—Decrease $0.5\%$	4.4979	0.4860	0.0473	2.0051
1c—Decrease $1.0\%$	4.4454	0.4920	0.0484	3.2647
2a—Increase 0.1%	N/A	N/A	N/A	N/A
2b—Increase 0.5%	4.6501	0.4686	0.0441	−1.6459
2c—Increase 1.0%	4.7184	0.4608	0.0427	−3.2844
3a—Bridge 10%	4.5820	0.4763	0.0455	−0.0123
3b—Bridge 15%	4.5839	0.4761	0.0454	−0.0578
3c—Bridge 20%	4.5892	0.4755	0.0453	−0.1850
4a—Cap 10%	4.5802	0.4766	0.0455	0.0309
4b—Cap 15%	4.5790	0.4767	0.0455	0.0597
4c—Cap 20%	4.5737	0.4773	0.0457	0.1869
True	4.5815	0.4764	0.0455

Strategy 2a yielded a small porosity error but a significantly finer mesh was required for numerical stability than that required for the other strategies. Results of this strategy were not discussed and denoted as not available (N/A).

summarizes the total pebble volume, calculated porosity, pore-scale hydraulic diameter, and the porosity error for each contact strategy relative to the true porosity. Some of the tested contact strategies yielded a significant porosity error relative to the true porosity. Decreasing or increasing the pebble diameter by $1\%$ yielded a porosity error of $3.2\%$. The bridging method with the smallest bridge diameter, equal to $10\%$ of the pebble diameter, showed the smallest porosity error of $0.012\%$. This strategy (Strategy 3a (bridge 10%)) was used as the baseline numerical result against which all the other strategies were compared. Strategy 2a, increasing the pebble diameter by 0.1%, yielded a small porosity error, but the simulations were not numerically stable unless a significantly finer mesh than required for the other strategies was used. Therefore, the results of this strategy were not discussed and denoted as not available where applicable.

STAR-CCM+ v16.02.009, a commercial finite-volume CFD software package developed by Siemens (Munich, Germany), was selected for this study. A steady Reynolds-averaged Navier–Stokes (RANS)-based approach to turbulence modeling was chosen for this effort given the size of the domain of interest and the available computational resources. This approach was chosen to maximize the accuracy of the numerical solution for conjugate heat transfer while balancing the relative computational costs to the porous media, large-eddy simulation, and direct numerical simulation methods.

2. Materials and Methods

2.1. Numerical Domain and Material Properties

Figure 1 presents the computational domain used in this study. A total of 100 pebbles were explicitly modeled using a $4\times 5\times 5$ structured array. Each pebble consisted of a fuel sphere surrounded by a fuel-free spherical shell. The fuel region represented a homogenization of TRISO particles within a graphite matrix. This region was defined as the power-producing volume with spatially constant volumetric heat sources (Figure 1b). Before making geometric modifications, the fuel region diameter, $d_{fuel}$, was 0.04 m, and the fuel-free region outer diameter, $d_p$, was 0.05 m. The heat transfer fluid inlet was located at the top of the domain, and the fluid outlet was located at the bottom. Three vertical planar slices were created to the domain to establish surfaces for symmetric boundary conditions. The resulting bed contained 45 whole pebbles and 55 sliced pebbles. The fourth vertical surface of the domain represented a radial reflector wall that was in contact with the pebbles. This surface was specified with a no-slip hydraulic boundary, a convective heat transfer coefficient of 10 W/(m$·$K), and a 250 °C ambient temperature to approximate the thermal pathway through a radial reflector and reactor vessel to a prototypical gas-cooled reactor cavity.

All boundary and interfacial conditions are specified in

Table 2. Boundary and interfacial conditions.

Region	Boundary	Condition	Value
Pebble fuel region	symmetry	symmetry plane	N/A
	fuel to fuel-free	contact interface	conformal conduction
Pebble fuel-free region	symmetry	symmetry plane	N/A
	fuel-free to fuel	contact interface	conformal conduction
	fuel-free to fluid	contact interface	conformal convection
	fuel-free to reflector	temperature	250 °C
Heat transport fluid	symmetry	symmetry plane	N/A
	inlet	inlet velocity 1	−3 ms−1
		inlet temperature	250 °C
		turbulence intensity	0.01
		turbulent viscosity ratio	10
	outlet	outlet pressure	$6\times {10}^6$ Pa
	fluid to reflector	temperature	250 °C
	fluid to fuel-free	contact interface	conformal convection

Adjusted for Case 2a (0 ms⁻¹), Case 2b (0 ms⁻¹), and Case 2c (−0.21 ms⁻¹).

for each region. The inlet turbulence parameters of turbulent kinetic energy, turbulent dissipation rate, turbulent viscosity, and turbulent thermal diffusivity were derived from the inlet specifications of 0.01 and 10 for turbulence intensity and turbulent viscosity ratio, respectively.

The material properties, except the fluid density, were specified at 6 MPa and 250 °C. These are independent of temperature and position within the computational domain and are presented in

Table 3. Specified temperature- and fluence-independent material properties.

Material Property	Pebble Fuel Region	Pebble Fuel-Free Region	Heat Transport Fluid
Density (kgm−3)	1700	1700	Ideal gas
Specific heat (Jkg−1K−1)	700	700	5200
Thermal conductivity (Wm−1K−1)	25	25	0.155
Dynamic viscosity ($\mathrm{Pa}\mathrm{s}$)	N/A	N/A	$1.98\times {10}^{-5}$
Molecular weight (kgkmol−1)	N/A	N/A	4.00
Turbulent Prandtl number [12]	N/A	N/A	0.7

. The fluid density was determined using a compressible flow solver and the ideal gas approximation. Also, the thermal contact resistance between pebbles was ignored, and the pebble chamfer bridges were treated as if they were made of the same material as the pebble fuel-free region. It is apparent and expected that pebble temperatures are sensitive to how the near-wall region of the fluid domain is treated. Conjugate heat transfer problems require both the velocity and thermal boundary layers to be appropriately resolved. The ratio of these boundary layer thicknesses can be estimated by the molecular Prandtl number, $P{r}_m$ represented as

$$P{r}_m={\displaystyle \frac{c_{p,f}\mu }{k_f}}={\displaystyle \frac{\mathrm{viscous}\mathrm{diffusion}\mathrm{rate}}{\mathrm{thermal}\mathrm{diffusion}\mathrm{rate}}}$$

(5)

where $c_{p,f}$ is the fluid-specific heat capacity, $\mu$ is the fluid-dynamic viscosity, and $k_f$ is the fluid thermal conductivity and yields a $P{r}_m$ of 0.664. This implies that the velocity boundary layer thickness is slightly less than the thermal boundary layer thickness. From a domain discretization standpoint, the conclusion can be made that the grid resolution requirements to accurately capture the near-wall velocity and temperature gradients are relatively similar.

2.2. Numerical Solvers

The heat-transfer and fluid-flow physics were coupled in a segregated manner of energy and flow transport. A detailed sensitivity study of turbulence models has been investigated but is not presented here in detail because discussions regarding RANS [13], unsteady Reynolds-averaged Navier–Stokes (URANS) [14], hybrid RANS/LES [15], and DNS [16] turbulence modeling comparisons already exist in the literature for body-centered cubic, face-centered cubic, structured, and unstructured. Most of the non-linear numerical methods generally predict similar mean velocity and temperature fields but differ significantly on the prediction of the root mean square (RMS) velocity and temperature fields with the bed. It has been suggested that the reason for poor RMS agreement is related to a significant amount of transitional flow affected by the pebble curvature.

However, validation data are lacking that are available to perform a thorough evaluation of turbulence modeling for the pebble-bed reactor application. Singular experimental configurations of heated beds that enable simultaneous measurement of temperature and time-resolved velocity fields are not easily accomplished. Prior isothermal efforts focused on particle image velocimetry (PIV) for time-resolved velocimetry [17,18] but lack temperature gradients.

In this study, RANS standard k-$ϵ$, realizable k-$ϵ$, and k-$\omega$ SST models were investigated. Relative differences in the predicted mean velocity fields between these three models could be decreased by

• Enabling the standard k-$ϵ$ cubic constitutive relationship to improve anisotopic turbulence prediction of secondary flows and streamline curvature [19].
• Enabing the k-$\omega$ SST cubic constitutive relationship [20,21], which adds non-linear functional representations of the strain and vorticity tensors.

Across these closure models, the volume-averaged mean velocity field and bed pressure drop spanned ±4% and the volume-averaged mean pebble fuel temperature spanned ±2%. The steady RANS realizable k-$ϵ$ turbulence model [22] with a second-order convection scheme, secondary gradients, and Wolfstein shear-driven two-layer all $y^{+}$ wall treatment [23] were selected. This decision was based on its wide applicability and the performance observed by Shams et al. [14] with acceptable agreement of mean fields with higher-order methods.

It is important to note that, regardless of the span of absolute values of velocity, pressure, and temperature fields, the conclusions regarding the geometric approach and selection of point-contact strategies remain unaffected.

The turbulent fluid motion and heat transfer are governed by mass continuity, solid energy conservation, fluid energy conservation, Navier–Stokes, and Poisson equations, provided as

$$\nabla ·\mathit{U}=0,$$

(6)

$${\displaystyle \frac{d}{dt}}{\int \int \int }_V{\rho }_s{c}_{p,s}TdV+{\oint }_A\rho c_{p,s}T{\mathit{v}}_{\mathit{s}}·d\mathit{a}=-{\oint }_A\dot{{\mathit{q}}_{\mathit{s}}^{″}}·d\mathit{a}+{\int \int \int }_V{S}_{u}dV$$

(7)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}{\int \int \int }_V{\rho }_{f}EdV+{\oint }_A{\rho }_{f}H\mathit{v}·d\mathit{a}& =\hfill \\ \hfill & -{\oint }_A\dot{{\mathit{q}}_{\mathit{f}}^{″}}·d\mathit{a}+{\oint }_A{\mathit{T}}_{\mathit{f}}·\mathit{U}d\mathit{a}\hfill \\ \hfill & +{\int \int \int }_V{\mathit{f}}_{\mathit{b}}·dV+{\oint }_A\sum _i{h}_i{J}_{i}d\mathit{a}\hfill \\ \hfill & +{\int \int \int }_V{S}_{u,f}dV\hfill \end{array}\end{array}$$

(8)

$${\displaystyle \frac{\partial U_i}{\partial t}}+U_j{\displaystyle \frac{\partial U_i}{\partial x_j}}=-{\displaystyle \frac{\partial P}{\partial x_i}}+\nu {\displaystyle \frac{{\partial }^2{U}_i}{\partial x_j\partial x_j}},$$

(9)

$${\displaystyle \frac{{\partial }^{2}p}{\partial x_i^2}}=-{\displaystyle \frac{\partial U_i{U}_j}{\partial x_i\partial x_j}}.$$

(10)

where $\mathit{U}$ is the fluid velocity vector, ${\rho }_s$ is the solid density, $c_{p,s}$ is the solid specific heat capacity, T is the solid temperature, $\dot{{\mathit{q}}_{\mathit{s}}^{″}}$ is the solid heat flux vector, ${\mathit{v}}_{\mathit{s}}$ is the solid convective velocity, $S_u$ is the solid volumetric heat source, ${\rho }_f$ is the fluid density, E is the fluid total energy, H is the fluid total enthalpy, $\dot{{\mathit{q}}_{\mathit{f}}^{″}}$ is the fluid heat flux vector, ${\mathit{T}}_{\mathit{f}}$ isthe viscous stress tensor, ${\mathit{f}}_{\mathit{b}}$ is the fluid body force vector, $h_i$ is the fluid enthalpy of component i, $J_i$ is the diffusive flux of component i, $S_{u,f}$ is the fluid volumetric heat source, P is the fluid internal pressure, t is time, $x_i$ is the spatial coordinate vectors of component i, and $\nu$ is the fluid kinematic viscosity.

The steady RANS-based approach to turbulence modeling was chosen for this study due to its relatively low computational cost to porous media when compared to that of LES and DNS. The temporal averaging methodology, or Reynolds decomposition process [24], suggested by Osborne Reynolds was to separate fluid variables into their time-averaged mean and fluctuating components. Applying this process to Equations (6)–(10) yields

$${\displaystyle \frac{\partial 〈U_i〉}{\partial x_i}}=0,$$

(11)

$${\displaystyle \frac{\partial 〈U_i〉}{\partial t}}+〈U_j〉{\displaystyle \frac{\partial 〈U_i〉}{\partial x_j}}+{\displaystyle \frac{\partial 〈u_i{u}_j〉}{\partial x_j}}=-{\displaystyle \frac{\partial 〈P〉}{\partial x_i}}+{\nu }_t{\displaystyle \frac{{\partial }^2〈U_i〉}{\partial x_j^2}},$$

(12)

$${\displaystyle \frac{{\partial }^2〈P〉}{\partial x_i^2}}=-{\displaystyle \frac{\partial 〈U_i〉}{\partial x_j}}{\displaystyle \frac{\partial 〈U_j〉}{\partial x_i}}-{\displaystyle \frac{{\partial }^2〈u_i{u}_j〉}{\partial x_i\partial x_j}}.$$

(13)

After the Reynolds decomposition of the Navier–Stokes equations, a closure problem is required to model the Reynolds stress tensor, $〈u_i{u}_j〉$. In this study, the applied closure method is the realizable $k-ϵ$ model containing a new transport equation for the turbulent dissipation rate $ϵ$ [22]. The transport equations for the kinetic energy k and the turbulent dissipation rate $ϵ$ are

$${\displaystyle \frac{\partial }{\partial t}}\left(\rho k\right)+\nabla ·\left(\rho k\overline{\mathit{v}}\right)=\nabla ·\left[\left(\mu +{\displaystyle \frac{{\mu }_t}{{\sigma }_k}}\right)\nabla k\right]+P_k-\rho (ϵ-{ϵ}_0)+S_k,$$

(14)

$${\displaystyle \frac{\partial }{\partial t}}\left(\rho ϵ\right)+\nabla ·\left(\rho ϵ\overline{\mathit{v}}\right)=\nabla ·\left[\left(\mu +{\displaystyle \frac{{\mu }_t}{{\sigma }_k}}\right)\nabla ϵ\right]+{\displaystyle \frac{1}{T_e}}{C}_{ϵ1}{P}_{ϵ}-C_{ϵ2}{f}_2\rho ({\displaystyle \frac{ϵ}{T_e}}-{\displaystyle \frac{{ϵ}_0}{T_0}})+S_e,$$

(15)

where $\overline{\mathit{v}}$ is the mean fluid velocity vector and ${\mu }_{t,f}$ is the fluid turbulent kinematic viscosity. ${\sigma }_k$, ${\sigma }_{ϵ}$, $C_{ϵ1}$, and $C_{ϵ2}$ are model coefficients, $P_k$ and $P_{ϵ}$ are production terms, $f_2$ is a damping function, $T_e$ is the large-eddy time scale, $T_0$ is the ambient source term time scale, $S_k$ and $S_{ϵ}$ are the user-specified source terms, and ${ϵ}_0$ is the ambient turbulence value in the source terms that counteracts turbulence decay.

For conjugate heat transfer of conduction and convection, the numerical method within STAR-CCM+ for adjacent fluid/solid cells create a fluid/solid contact interface and conserves the total heat flux. The heat flux is linearized with the following form:

$$q_0^{″}=A_0+B_0{T}_{f0}+C_0{T}_{w0}+D_0{T}_{w0}^4$$

(16)

where A, B, C, and D are net wall heat flux coefficients, $T_{f0}$ is the bulk fluid temperature, and $T_{w0}$ is the interfacial temperature. The heat transfer coefficient is defined locally based on the local flow conditions (fluid density, specific heat capacity, wall shear stress, and local Reynolds number), local fluid temperature, and local solid temperature. And this results in a local heat transfer coefficient for each cell face along the interface [23].

In Case 2, PLOFC, thermal radiative heat transfer was also modeled using the surface radiation exchange model (S2S) [23]. The heat transport fluid was assumed to be non-participating with no internal absorption, emission, or scatter (100% transparent to thermal radiation). The pebble wetted surfaces were assigned as the participating faces for thermal radiation, and these were used to calculate the appropriate view factors. No radiative heat transfer was modeled within the solid pebbles (100% opaque).

The radiation exchange between surfaces can be accounted for by tracking the amount of energy emitted from and absorbed from pairs of surface mesh faces as

$$P_{1-2}=i_1^{\prime }d{S}_{1}cos\left({\beta }_1\right)({\displaystyle \frac{d{S}_{2}cos\left({\beta }_2\right)}{L^2}})$$

(17)

where $P_{1-2}$ is the power emitted from face 1 and absorbed by face 2, $i_1^{\prime }$ is radiation intensity flux per unit area and per unit steradian, $\beta$ is the angle between faces, and L is the distance between faces.

Conservation of radiative energy is balanced only when using a closed set of surfaces. The total energy transfer can be calculated by summing over all faces that represent the closed set. For face i, this becomes

$$Q_{I,i,\lambda }=S_i{I}_{i,\lambda }=\sum _{j=1}^{N_p}{\displaystyle \frac{F_{j-i,\lambda }}{(1-{\tau }_{i,\lambda }-{\rho }_{s,i,\lambda })}}{J}_{d,j,\lambda }{S}_j+{\displaystyle \frac{F_{e-i,\lambda }}{(1-{\tau }_{i,\lambda }-{\rho }_{s,i,\lambda })}}{J}_{e,\lambda }{S}_e$$

(18)

where $S_i$ and $S_j$ are the indexed face surface areas, $I_{i,\lambda }$ is the indexed radiation flux, $N_p$ is the total number of faces, $F_{j-i,\lambda }$ is the view factor from face j to i, $J_{d,j,\lambda }$ is the diffuse component of radiation leaving face j, $F_{e-i,\lambda }$ is the view factor from the environment to face i, ${\tau }_{i,\lambda }$ is the transmissivity, and ${\rho }_{s,i,\lambda }$ is the specular component of radiation. Other radiative assumptions included

• Each pebble surface mesh face was used in the view factor calculation;
• Gray radiation with no wavelength dependence;
• Diffuse radiation with no angular emission dependence;
• The effective radiation temperature of the environment was specified as 250 °C;
• A pebble surface emissivity of 0.8; and
• Satisfaction of Kirchhoff’s Law, such that emissivity + reflectivity + transmissivity = 1. This forces the surface reflectivity to 0.2 because the transmissivity is 0.

2.3. Mesh Sensitivity

Tetrahedral, hexahedral, and polyhedral mesh generation approaches were tested for each contact strategy. The bridging and diametral increase strategies resulted in a singular, complex body of all the fuel-free spherical shell regions. The parallel load distribution step required more computational time to decompose this body across the compute nodes. The capping and diametral decrease strategies resulted in many simple bodies. These geometries were significantly faster for interfacial definition (imprinting) and parallel discretization.

Fully conformal meshes were generated using polyhedral and prism layer meshing operations. The Surface Remesher performed surface vertex re-tessellation of the 3D CAD to optimize surface faces based on the target edge length and proximity refinements. In the Polyhedral Mesher process, a tetrahedral mesh is first generated for the input surface. Second, the polyhedral mesh is generated from the underlying tetrahedral mesh. The Prism Layer Mesher generated a subsurface to extrude a set of prismatic cells from the solid–fluid interfaces into the core mesh of the fluid region.

The sensitivity of the discretization on fluid flow and heat transfer behavior was quantified by performing mesh sensitivity studies for Cases 1–2. The Case 1 study is discussed in detail for each contact strategy. The diametral increase and bridging strategies had similar grid requirements. The diametral decrease and capping strategies required additional grid refinement to improve the prediction of the local velocity gradient in the small fluid pathways near the contact point. This refinement was also needed to achieve mesh independence for the pressure drop across the bed. For brevity, only the results from the bridging strategy are presented here. Coarse, medium, and fine meshes were generated using base cell sizes of 10 mm, 2 mm, and $0.7$ mm, respectively. The resulting numbers of cells for each mesh were $2.4\times {10}^6$, $7.4\times {10}^6$, and $41.0\times {10}^6$. Figure 2 presents cross-sections of the bed for each mesh resolution.

To evaluate the numerical discretization error, estimate the rate of grid convergence, and determine an acceptable mesh density for this problem, the widely accepted Richardson extrapolation technique [25] was used. Quantities of interest (QOIs) were selected to quantify the grid convergence:

• Volume-averaged pebble fuel temperature, $T_{avg,fuel}$;
• Maximum pebble fuel temperature, $T_{max,fuel}$;
• Volume-averaged fluid temperature, $T_{avg,fluid}$;
• Volume-averaged fluid velocity magnitude, $v_{avg,fluid}$;
• Maximum fluid velocity magnitude, $v_{max,fluid}$;
• Surface-averaged fluid temperature on plane Z2 (z = 0.1 m), $T_{avg,fluid,trans}$;
• Surface-averaged fluid velocity magnitude on plane Z2 (z = 0.1 m), $v_{avg,fluid,trans}$;
• The pressure drop across the entire bed, $\Delta P_{bed}$.

The convergence of each QOI, ${\varphi }_i$, is shown in Figure 3 as a function of a 3D unstructured grid representative cell size, h, and normalized by the extrapolated QOI with an infinite grid density, ${\varphi }_{extrap}$. For the coarse, medium, and fine meshes, $h_i$ was defined as 2.15, 1.50, and 0.84 mm, respectively, using the definition of

$$h={\left[{\displaystyle \frac{1}{N}}\sum _{i=1}^N(\Delta V_i)\right]}^{1/3}.$$

(19)

Macroscopic QOIs, such as the volume-averaged metrics, had similar behavior and rates of convergence across the various contact strategies. For the coarse mesh, all relevant QOIs were within 3% of an extrapolated solution with an infinite grid density, except for the velocity magnitude selected at an arbitrary point in the bed. This discretization error was reduced to nearly 1% for all QOIs for the fine mesh.

It is important to note that $v_{max,fluid}$, a point-wise QOI, was selected at an arbitrary coordinate in the bed and was not held spatially constant. Also, the polyhedral mesh refinement did not preserve cell centroids and thus the interpolation scheme and shifting cell centers induce additional uncertainty onto this point-wise term. The Richardson extrapolation technique returned a higher-than-expected order of convergence for this QOI, which highlights the sensitivity of the method and emphasizes the care one must take in selecting QOIs. The bed pressure drop, $\Delta P_{bed}$, converged oscillatorily for Strategies 2 and 3. This may be due to the unsteady nature of the solution of flow around spherical bodies in combination of flow around cylindrical and near-cylindrical bodies.

Therefore, for $v_{max,fluid}$ and $\Delta P_{bed}$, the Richardson extrapolation was not applied. Instead, the average and maximum difference between the solutions were calculated and the difference multiplied by 3 was deemed as the grid convergence uncertainty. For $v_{max,fluid}$, this resulted in an uncertainty of ±7%. For $\Delta P_{bed}$, this resulted in an uncertainty of ±3%.

Finally, the medium mesh was selected as the working grid size for further study because most QOIs were within 2% of the extrapolated solution while requiring only 18% of the fine mesh volumetric element count.

3. Results

3.1. Case 1: Full Power

In Case 1, each contact strategy was investigated regarding the full-power condition.

Table 4. Relative differences in QOIs computed from various pebble contact strategies for Case 1 (%).

Strategies	Porosity Error [%]	ΔP	${\mathit{v}}_{\mathit{avg},\mathit{fluid}}$	${\mathit{T}}_{\mathit{max},\mathit{fuel}}$	${\mathit{T}}_{\mathit{avg},\mathit{fuel}}$	${\mathit{T}}_{\mathit{avg},\mathit{fluid}}$
1a—Decrease 0.1%	0.3284	−0.3	−0.4	−0.1	0.0	0.0
1b—Decrease 0.5%	2.0051	−22.4	−4.3	0.2	0.2	0.0
1c—Decrease 1%	3.2647	−33.0	−7.3	0.3	−0.1	0.0
2a—Increase 0.1%	N/A	N/A	N/A	N/A	N/A	N/A
2b—Increase 0.5%	−1.6459	12.4	2.8	0.2	0.0	0.0
2c—Increase 1%	−3.2844	21.1	5.1	0.2	0.0	0.0
3a—Bridge 10%	−0.0123	0.0	0.0	0.0	0.0	0.0
3b—Bridge 15%	−0.0578	1.8	−0.2	0.2	0.1	0.0
3c—Bridge 20%	−0.1850	3.1	0.5	0.6	0.3	0.0
4a—Cap 10%	0.0309	−13.1	−1.9	0.0	0.0	−0.1
4b—Cap 15%	0.0597	−22.8	−3.6	0.2	−0.4	0.0
4c—Cap 20%	0.1869	−27.3	−4.4	−0.4	−0.3	0.0
GCI (%)—Bridge 10%	N/A	3	0.1	3.4	0.6	0.1

presents the QOIs computed for each contact strategy, including the following: the pressure drop across the entire bed, $\Delta P_{bed}$, the maximum velocity magnitude within the bed, $v_{max,fluid}$, the volume-averaged velocity magnitude within the bed, $v_{avg,fluid}$, the maximum pebble temperature, $T_{max,fuel}$, the volume-averaged pebble temperature, $T_{avg,fuel}$, and the volume-averaged fluid temperature within the bed, $T_{avg,fluid}$.

Strategy 3a (bridging 10%) was used as the baseline reference against which all the other contact strategies were compared because it had the smallest true porosity error of −0.01%. Figure 4 presents the pebble temperature, fluid temperature, and fluid velocity fields for the baseline Strategy 3a.

In general, the relative agreement of the other contact strategies decreased as the porosity error increased. The hydraulic QOIs were more sensitive than the thermal QOIs. The sensitivity of the contact strategy on the hydraulic QOIs was approximately 10 times greater than the observed differences when compared to the ${\varphi }_{extrap}$ in the mesh convergence study. The sensitivity of the contact strategy on the thermal QOIs was approximately equal to the observed differences when compared to the ${\varphi }_{extrap}$ in the mesh convergence study.

In Strategy 1, reducing the pebble volume, the Strategy 1a approach of a 0.1% diametral reduction performed the best: less than $1\%$ differences were observed for all the QOIs when compared to the baseline. On the other hand, Strategy 1c, with the error of porosity at $3.264\%$, had the highest difference in QOIs; it significantly under-predicted the pressure drop and average velocity magnitude at $33.0\%$ and $7.3\%$, respectively. For Strategy 2, this approach of increasing the pebble diameter was found to overestimate the pressure drop and average velocity magnitude. Only increasing the diameter by 0.1% yielded a small porosity error, but the simulations were not numerically stable unless a significantly finer mesh than required for the other strategies was used. A $1\%$ diametral increase yielded a porosity error of $-3.28\%$ and large relative differences in pressure drop ($21.1\%$) and average velocity magnitude ($5.1\%$). With respect to Strategy 3, changing the bridge diameter from $10\%$ to $20\%$ of the pebble diameter had a small impact. The largest differences were found to be less than $3.2\%$. Strategy 4 underestimated the pressure drop and slightly underestimated the average fluid velocity, but it still yielded very similar pebble and fluid temperatures compared to the baseline.

Figure 5 illustrates the surface temperature contours of the numbered pebbles along the views from a side view ($xz$-plane), a top view ($xy$-plane), and a bottom view ($xy$-plane) that were obtained using Strategies 1a, 2b, 3a, and 4a. Recall that Strategies 1 and 2 involve a simple modification of the pebble diameter, which is easy to implement for very large beds. On the other hand, Strategies 3 and 4 require considerable effort to produce the bridging between the contacting pebbles but were found to have lower porosity errors. Slight temperature differences are visible on the pebble downstream halves when viewed from the $xy$-plane bottom view.

Figure 6 depicts line integral convolution streamlines and velocity magnitudes along the $xz$- and $xy$-planes for Strategies 1a, 2b, 3a, and 4a. Among these contact strategies, the peaks of the velocity magnitude obtained from Strategies 1a and 3a appear the most similar for both the $xz$- and $xy$-planes. When comparing the velocity contours along the $xz$-plane (side view), the spatial positions of the velocity maxima were found to switch sides based on the contact strategy. For example, see pebble group 4-5-10-11. For the subchannels located in the vicinity of the reflector wall—that is, those that were created by pebbles 10-4-13 in the $xz$-plane and by pebbles 1-4-7 in $xy$-plane—all the tested pebble contact strategies revealed similar flow patterns. For instance, in each subchannel of pebbles 1-4-7, two pairs of counter-rotating vortices impinge on the pebble surfaces in each strategy.

3.2. Case 2: PLOFC

In Case 2, the four contact strategies were again investigated and roughly correspond to the PLOFC with a reactor trip event. This case was motivated by an objective to better understand the contact strategy under the separate effects of each heat transfer mode of conduction, convection, and radiation. For this case, the pebbles were assumed to have a temporally constant power density equal to 7% of the full-power density, which corresponds to approximately 20 s after the reactor trip [10].

Table 5. Relative differences in QOIs computed from each pebble contact strategy for Cases 2a and 2b (%).

Strategy	Case 2a			Case 2b
	${\mathit{T}}_{\mathit{max},\mathit{fuel}}$	${\mathit{T}}_{\mathit{avg},\mathit{fuel}}$	${\mathit{T}}_{\mathit{avg},\mathit{fluid}}$	${\mathit{T}}_{\mathit{max},\mathit{fuel}}$	${\mathit{T}}_{\mathit{avg},\mathit{fuel}}$	${\mathit{T}}_{\mathit{avg},\mathit{fluid}}$
1a—Decrease $0.1\%$	13.8	6.5	5.8	1.4	0.7	0.2
1b—Decrease $0.5\%$	20.2	9.7	8.9	1.6	0.7	0.4
1c—Decrease $1.0\%$	25	12.1	11.5	1.5	0.7	0.6
2a—Increase 0.1%	N/A	N/A	N/A	N/A	N/A	N/A
2b—Increase 0.5%	−1.5	−0.9	−1.2	0	0	−0.2
2c—Increase 1.0%	−4.3	−2.5	−3.1	−0.6	−0.2	−0.4
3a—Bridge 10%	0.0	0.0	0.0	0.0	0.0	0.0
3b—Bridge 15%	−2.9	−1.4	−1.2	−0.7	−0.3	−0.1
3c—Bridge 20%	−4.7	−2.3	−2	−1.3	−0.6	−0.2
4a—Cap 10%	15.6	7.8	8.1	0.5	−0.5	−1.3
4b—Cap 15%	19.6	9.6	9.7	0.6	−0.5	−1.3
4c—Cap 20%	22.8	11.1	11.1	0.8	−0.4	−1.2

and

Table 6. Relative differences in QOIs computed from various pebble contact strategies for Case 2c (%).

Strategies	$\Delta \mathit{P}$	${\mathit{v}}_{\mathit{avg},\mathit{fluid}}$	${\mathit{T}}_{\mathit{max},\mathit{fuel}}$	${\mathit{T}}_{\mathit{avg},\mathit{fuel}}$	${\mathit{T}}_{\mathit{avg},\mathit{fluid}}$
1a—Decrease 0.1%	−1.2	−0.4	0.0	0.0	0.0
1b—Decrease 0.5%	−8.6	−2.3	0.1	0.0	0.0
1c—Decrease 1%	−19.1	−5.1	0.2	0.1	0.0
2a—Increase 0.1%	N/A	N/A	N/A	N/A	N/A
2b—Increase 0.5%	9.8	2.3	−0.1	−0.1	0.0
2c—Increase 1%	18.4	4.4	−0.1	−0.2	0.0
3a—Bridge 10%	0.0	0.0	0.0	0.0	0.0
3b—Bridge 15%	−0.4	0.2	0.0	0.0	0.0
3c—Bridge 20%	−0.9	0.5	0.0	0.1	0.0
4a—Cap 10%	−0.5	0.0	0.3	0.2	0.0
4b—Cap 15%	−2.3	−0.5	0.3	0.1	0.0
4c—Cap 20%	−7.2	−3.4	0.6	−0.1	0.0

summarize the relative differences in the QOIs for each subcase. For Case 2a (conduction), the results are naturally divided into two groups based on whether or not the pebble bed was conductively connected (e.g., Strategies 2 and 3) or conductively isolated (e.g., Strategies 1 and 4). The conductively connected results had less than a 5% difference in the fuel and fluid temperatures, whereas the conductively isolated results differed significantly and diverged as the pebble volumes continued to decrease. The most different results were yielded by Strategies 1c and 4c, at approximately 25% and 22.9% differences, respectively, when predicting the maximum fuel temperature compared to the baseline results of Strategy 3a. As a consequence, the maximum fuel temperatures predicted by Strategies 1 and 4 were found to be considerably higher than those obtained from Strategies 2 and 3, where the pebble contacts are maintained. In addition, the fluid temperatures estimated using Strategies 1 and 4 were observed to be greater than those computed using Strategies 2 and 3; up to a 12.3% difference in $T_{avg,fluid}$ was predicted using Strategy 1c.

For Case 2b with conduction and radiation heat transfer modes, the temperature differences decreased significantly in both the solid and fluid regions. Strategy 1 slightly over-predicted the maximum fuel temperature within a 1.6% difference and fluid temperature with less than a 0.7% difference compared to the baseline results. The results of Strategy 2 were similar to the baseline, with differences of less than 0.7%. For Strategy 3, increasing the bridge size yielded greater differences in both the fuel and fluid temperatures compared to the baseline bridge 10% d_p method. The temperatures computed using Strategy 4 varied from 0.4% to 0.8%.

Figure 7 presents the Case 2c pebble surface temperatures along the views from the $xz$-plane (side view), $xy$-plane (top view), and $xy$-plane (bottom view) using Strategies 1a, 2b, 3a, and 4a. The pebble surface temperature contours from Strategy 2b were found to be lower than those from the baseline and other strategies. Table 6 highlights that the fluid velocity is the highest for Strategy 2, which explains the lower surface temperature. Nevertheless, the maximum and volume-averaged temperatures are nearly identical across the strategies.

Although the temperatures were not sensitive to the contact strategy, much larger variances were observed in the predictions of the maximum velocity—nearly 7.5% using Strategy 4c—and volume-averaged velocity, at 5.1% using Strategy 1c, within the bed. It can be concluded that any modification to the bed that increases or decreases the porosity, even slightly, results in larger differences in fluid velocity. This trend was also observed in the predictions of pressure drops from all the tested strategies. For instance, Strategy 1a (decrease 0.1%) with a 0.33% porosity error and Strategy 3a (decrease 1%) with a 3.26% porosity error yielded differences in pressure drop predictions of 1.2% and 19.1%, respectively, compared to the baseline. Similarly, varying the cap sizes from 10% (0.031% porosity error) to 20% (0.187% porosity error) caused the differences in the pressure drops to increase from 0.52% to 7.2%, respectively.

4. Discussion

This study numerically investigated the effects of various contact strategies on the flow and heat transfer behavior within a structured bed of 100 explicitly modeled pebbles. The pebbles were modeled as perfect spheres, and each had an inner spherical fuel region and an outer spherical shell fuel-free region. Four contact strategies and two thermal hydraulic conditions were considered. The strategies to avoid point-contact singularities included decreasing and increasing the pebble diameter, along with capping and bridging the pebble surfaces near the contact region. The two thermal hydraulic conditions were full-power operation (Case 1) and pressurized loss of forced cooling (Case 2).

The most desirable contact strategies are those that prioritize minimizing the porosity error as much as possible while providing a thermal connection to enable specifying an effective point-contact thermal resistance. Therefore, the recommended contact strategies and sizing selection are expected to hold for unstructured beds. Also, a study investigating externally heated unstructured beds reached similar conclusions as this effort with internally heated pebbles regarding the opportunity cost of each strategy and sizing recommendations [1].

Strategy 3a, bridging with a cylinder equal to 10% of the pebble diameter, was selected as the baseline strategy because it addressed the contact singularity while minimizing the geometric changes that affect the bed porosity.

For Case 1 (full power), the hydraulic behavior was sensitive to the various contact strategies, with a 30% difference in the predicted pressure drop (3% GCI) and a 4% difference for the volume-averaged velocity magnitude (0.1% GCI). The thermal behaviors within the pebbles and fluid were not sensitive to these changes as a difference of less than 0.5% (0.6% GCI) was observed across the strategies.

The strategy most similar to the baseline (Strategy 3a, bridge 10%) was Strategy 1a, with a pebble diametral reduction of 0.1%. However, the small fluid volumes created from the diametral reduction required an additional +12% of the volumetric elements to adequately resolve the RANS flow field.

A benefit of Strategy 3 is the simplification of specifying the interfacial thermal contact resistance between pebbles. This is because the bridges can be kept as separate regions during the model-building process. However, this bridge building is not trivial for large and unstructured beds. The contact patch definitions required for conformal interfaces become increasingly difficult for large-bed model building with bridging. These can be avoided with diametral reduction (Strategy 1) and chamfering (Strategy 4).

For Case 2 (PLOFC), large differences between the contact strategies were observed in the predicted pebble temperature when only considering thermal conduction. The difference in the maximum pebble temperatures for various strategies was 23% in Case 2a. When including radiation in Case 2b, the difference decreased to 2%. When adding 7% flow in Case 2c, the thermal differences decreased further to 0.3% for the various strategies. The hydraulic behavior’s sensitivity to the contact strategies in Case 2c was similar to the sensitivity observed in Case 1, where Strategy 1a was in close agreement with the baseline.

If the objective is to accurately predict the pressure drop and velocity distribution, Strategies 1a, 3, and 4a are recommended. Any larger diametral changes or large bridge sizes have too large of an impact on the bed porosity (>2% porosity error for global strategies). To accurately predict the temperature distribution within the fuel pebbles, any of the studied point-contact strategies are suitable when neglecting the thermal contact resistance. Strategy 3 simplifies the implementation of the interfacial thermal contact resistance between pebbles. Nevertheless, the geometry modifications of Strategy 1a can be achieved significantly faster for large and unstructured pebble beds.

Future efforts will focus on generating the experimental data of the temperature, pressure, and velocity fields for each contact strategy. For validation purposes, the data will then be compared with the numerical results of this paper.