Improved Lumped Models for Transient Combined Convective and Radiative Cooling of Multilayer Spherical Media

In this work, we studied the transient combined convection and radiation of multilayer spherical media with volumetric heat generation, extending the previous work on the particular case of a spherical body subjected to radiative cooling. The proposed lumped models were obtained through two-point Hermite approximations for the average temperature and heat flux in each layer. For the average temperature, the plain trapezoidal rule (H0,0 approximation) was employed in all layers, except for the innermost layer, where the second-order two-side corrected trapezoidal rule (H2,1 approximation) was utilized. For the heat flux, the plain trapezoidal rule (H0,0 approximation) was employed for all the layers.The transient heat conduction in a TRISO-coated fuel particle being composed of five layers (namely, fuel kernel, buffer of porous carbon, inner pyrocarbon, silicon carbide, and outer pyrocarbon) was analyzed using the proposed lumpedmodels, the results of whichwere verified by comparisonwith the finite difference solution of the original distributed parameter model. Parametric studies were conducted to examine the effects of the dimensionless heat generation rate, the radiation-conduction parameter, and the Biot number on the temporal variations of the average temperatures.


Introduction
Due to additional advantage of combining the good mechanical properties and high thermal performance of different layers, multilayer structures have found widely applications in various technological areas, such as automotive, space, chemical, civil, and nuclear industries.Therefore, one of the considerations for the design of multilayer structures is to accurately and efficiently predict the temporal variation of the temperature of each layer.For instance, the thermal analysis of the triple-isotropic (TRISO) coated fuel particles consisting of fuel kernel (mainly UO 2 ), porous graphite buffer layer, inner dense pyrolytic carbon (IPyC) layer, SiC layer, and outer dense pyrolytic carbon (OPyC) layer is essential for safe and economic operation of high temperature gas-cooled reactors (HTGR) [1][2][3].
In various weight-saving and impact-property demanding applications, one often employs multilayer spherical composites as advanced structural components.Many analytical and numerical approaches have been proposed to obtain the transient thermal response of multilayer spherical media.Özisik [4] gave a detailed review of one-dimensional multilayer heat conduction problems, mainly focusing on orthogonal expansions by using Greens functions and Laplace transform techniques.By means of the transfer function method, Virseda and Pinazo [5] established a methodology to determine heat flux and inner/outer temperatures for spherical solid multilayer products.de Monte [6] derived a "natural" analytical approach for solving transient temperature response of multilayer composite spheres.Chen et al. [7] introduced the frequency-domain regression method to calculate the thermal response factors and the conduction transfer function coefficients of multilayer spherical structures.Ramadan [8] presented a semianalytical solution procedure for transient heat transfer in composite mediums consisting of multilayers within the framework of dual phase lag model, the solutions of which can be applicable to the classical Fourier heat diffusion, hyperbolic heat conduction, phononelectron interaction, and phonon scattering models.With the separation of variables and the recursive Thomas algorithm, Norouzi et al. [9] provided an exact analytical solution for steady conductive heat transfer in multilayer spherical fiber reinforced composite laminates.By introducing the Laplace transform and an approximate inverse Laplace transform, Lu and Viljanen [10] presented an analytical approach to heat conduction in a layered sphere subject to a time-dependent temperature change.Singh et al. [11] and Jain et al. [12] proposed a closed form analytical double-series solution for the multidimensional unsteady heat conduction problem in polar coordinates with multiple layers in the radial direction.As an extension of the previous work, Jain et al. [13] presented an analytical series solution for transient boundary-value problem of heat conduction in - spherical coordinates, which is applicable in spherical or part-spherical multilayer geometries such as spherical cone, hemisphere, spherical wedge, and full sphere.Singh et al. [14] presented an approach based on the finite integral transform method to solve the asymmetric heat conduction problem in a multilayer annulus with time-dependent boundary conditions and heat sources.
In this work, we present improved lumped models for transient combined convective and radiative cooling of multilayer spherical media, which is a more general case for heat conduction problem of spherical structures, extending the previous work on the particular case of a spherical body subjected to radiative cooling [16].The lumped model approach is useful to simplify otherwise complex differential heat equation, which treats each layer as a single lump with a single temperature.Although comprehensive computational techniques such as finite difference, finite volume, and finite element methods can be applied to solve distributed parameter formulations, the lumped model involving fewer state variables and consequently less equations to be solved is usually preferred by engineers to predict and control the system behavior.The proposed lumped models are obtained through two-point Hermite approximations for the average temperature and heat flux in each layer [17,18].This approach has been shown to be efficient in a great variety of multilayer applications [19][20][21][22][23][24].For the average temperature, the plain trapezoidal rule ( 0,0 approximation) is employed in all layers, except for the innermost layer, where the second-order two-side corrected trapezoidal rule ( 2,1 approximation) is utilized.For the heat flux, the plain trapezoidal rule ( 0,0 approximation) is employed for all the layers.The transient heat conduction in a TRISO-coated fuel particle is analyzed using the proposed lumped models, the results of which are verified by comparison with the finite difference solution of the original distributed parameter model.A parameter study is then performed to investigate the effects of the heat generation rate, the radiation-conduction parameter, and the Biot number on the time-temperature history of each layer in the TRISO particle.

Mathematical Formulation
Consider one-dimensional transient heat conduction in multilayer spherical composite shells consisting of  concentric layers in spherical coordinates, subjected to convective heat transfer at the inner side and combined convective and radiative heat transfer at the outer side, as shown in Figure 1.The layers are labeled as 1 to  from the innermost layer to the superficial layer.Let  be the coordinate perpendicular to the layers,   ,  = 1, 2, . . .,  represent its value at the inner surface of each layer, and  +1 be the coordinate value at the outer surface of th layer.It is assumed that the thermophysical properties of the layers are homogeneous, isotropic, and independent of the temperature.The volumetric rate of heat generation in th layer is   (, ).The adjacent layers are assumed to be in perfect thermal contact.Initially, all the layers are at a specified uniform temperature  0 .Convective heat transfer occurs at the inner boundary surface  =  1 and the outer boundary surface  =  +1 , with constant heat transfer coefficients ℎ 1 and ℎ +1 , to environmental fluids with constant temperatures  ∞ .Besides, a constant radiation sink temperature   is specified at the outer boundary surface.
The mathematical formulation of the one-dimensional transient heat conduction problem is given by with the following boundary and interface conditions: and the initial conditions for each layer where   (, ) is the temperature in the th layer,  is the time,   (=   /    ) is the thermal diffusivity,   is the thermal conductivity,   is the density,   is the specific heat,  is the surface emissivity, and  is the Stefan-Boltzmann constant.
It should be noted that in general the environmental fluid temperature  ∞ differs from the radiation sink temperature   .It is convenient to introduce the adiabatic surface temperature   , defined by The boundary condition ( 5) can be rewritten with use of the adiabatic surface temperature The mathematical formulation given by the system of ( 1)-( 4), (6), and ( 8) can now be expressed in dimensionless form as follows: where the dimensionless parameters are defined by Let   represent the dimensionless thickness for each layer,

Lumped Models
Let us introduce the spatially average dimensionless temperature of the th layer as follows: where the heat source term is defined as Equations ( 18) are an equivalent integrodifferential formulation of the original mathematical model, (9), with no approximation involved.In classical lumped-differential analysis, the boundary temperatures are assumed to be the same as the average temperatures.The approach is limited to low Biot numbers.In an attempt to enhance the approximation approach of the classical lumped model, we develop improved lumped models by providing better relations between the boundary temperatures and the average temperatures, based on Hermite-type approximations for integrals that define the average temperatures and the heat fluxes.The general Hermite approximation for an integral, based on the values of the integrand and its derivatives at the integration limits, is written in the following form [18]: where () and its derivatives  (]) () are defined for all  ∈ (, ).It is assumed that the numerical values of  (]) () for ] = 0, 1, . . .,  and  (]) () for ] = 0, 1, . . .,  are available.The general expression for the  , approximation is given by where ℎ =  −  and We first employ the plain trapezoidal rule to the integrals of both average temperatures for  = 2, 3, . . .,  and average heat fluxes for  = 1, 2, . . .,  ( 0,0 / 0,0 approximation), in the form Now, note that, for the average temperature of the innermost layer, if  0,0 or  1,0 approximation is employed, the average temperature is expressed by one term only ( 1 (, )| = 2 ), which will bring the calculation error.Therefore, applying  2,1 approximation to the integral of the definition of the average temperature of the innermost layer and supposing that we are studying a compact sphere ( 1 = 0), we have  10)-( 13) and ( 23)-( 24) by using symbolic computation software such as Mathematica and then used to close the ordinary differential equations (18) for the average temperatures  av ,  = 1, 2, . . ., , to be solved with the initial condition  av ()    =0 = 1,  = 1, 2, . . ., . (25)

Numerical Results and Discussions
As a typical thermal analysis problem of the multilayer spherical media, the temperature transients of the TRISOcoated fuel particles are often paid attention to due to their importance for design decision and safety assessment of the high temperature gas-cooled reactor.In this section, the verification of the proposed lumped models and the parametric study for the temperature transients of the TRISO particle are presented.

Verification of the Solution.
Transient heat conduction in 5-layer TRISO-coated gas-reactor particle fuel with combined convective and radiative cooling is analyzed to illustrate the applicability of the proposed lumped models.To demonstrate the validity and accuracy of the proposed approach, the solutions of the lumped model with Φ  = 0 and  rc = 0 are presented in graphical forms in comparison with a reference finite difference solution of the original distributed model, ( 9)-( 14).Based on the convergence analysis of the finite difference discretization, the initial boundary-value problem defined by ( 9)-( 14) is solved by using an implicit finite difference method, with 10-node mesh in spatial discretization for each layer and a dimensionless time step of 0.001.After the finite difference solution is obtained at each time step, the trapezoidal rule is used to integrate numerically the temperature in space to obtain the spatially average temperature.The thermal conductivities, the densities, the specific heat, and the geometries of the TRISO particle are shown in Table 1 [15], where  UO 2 is the diameter of the UO 2 kernel and  Buf ,  IPyC ,  SiC , and  OPyC are the thicknesses of the buffer, the IPyC, the SiC, and the OPyC layer, respectively.
The Biot number at the outer surface is specified to be 1.0.2.The initial boundary-value problem defined by ( 9)-( 14) is solved by using an implicit finite difference method.The good agreement can be found for the calculated averaged transient temperatures of each layer in the TRISO particle including UO 2 kernel, buffer layer, IPyC layer, SiC layer, and OPyC layer, as shown in Figure 2.

Parametric Study.
In this section, using the feasible approach, different values of the dimensionless heat generation rate Φ 1 , the radiation-conduction parameter  rc , and the Biot number Bi are chosen to assess their effects on the temperature history of each layer in the TRISO particle.
Other parameters including the thermal conductivities, the densities, the specific heat, the geometries, the environmental fluid temperature, the adiabatic surface temperature, and the initial temperature are the same as those employed in the previous section.lumped model are shown in Figure 3, where the dimensionless heat generation rate Φ 1 = 4.0, 8.0, 12.0, and 16.0, respectively.The radiation-conduction parameter  rc is equal to 5.0, and the Biot number Bi is specified to be 3.0.It can be observed that as the time increases, the transient average temperature for each layer in the TRISO particle decreases and approaches the steady-state conditions at around  = 0.2, which descends to 0.30, 0.40, 0.50, and 0.60, respectively.At a certain time, the temperature of the UO 2 kernel is always higher than the other four layers no matter the transient or steady-state phases in the cooling of the TRISO particle.In addition, the temperature in the UO 2 layer becomes higher as the heat generation rate grows for both the transient stage and the steady stage, while the temperatures in the other layers do not obviously change.

The Effect of the Radiation-Conduction
Parameter.The radiation-conduction parameter,  rc , is one of the governing parameters for convective cooling.The temporal variations of the average temperatures predicted by the lumped model are shown in Figure 4, where the dimensionless radiationconduction parameter  rc = 5.0, 10.0, 15.0, and 20.0, respectively.The heat generation rate Φ 1 is equal to 6.0, and the Biot number Bi is specified to be 3.0.For each case, the transient average temperature of all the layers decreases with the time.The higher  rc causes the transient average temperature to be lower, but the effect on  av of each layer is not quite remarkable as the one of the heat generation rate.

The Effect of the Biot Number.
The Biot number reflects the relation between convective and conductive flux.The temporal variations of the average temperatures predicted by the lumped model are shown in Figure 5, where the Biot number Bi = 3.0, 6.0, 9.0, and 12.0, respectively.The heat generation rate Φ 1 is equal to 6.0, and the dimensionless radiation-conduction parameter  rc is specified to be 8.0.It can be seen from the history of the transient average temperature calculated that increasing the value of Bi leads to a decrease in temperature of all the layers.

Conclusions
Based on two-point Hermite approximations for integrals, improved lumped parameter models are developed for the transient combined convective and radiative cooling of multilayer spherical media.Considering the TRISO particle as an example,  2,1 approximation is employed for the innermost layer, while  0,0 approximation is used for the other layers.The parametric studies show that as the time increases, the transient average temperature for each layer in the TRISO particle decreases and approaches the steady-state conditions.In addition, increasing the values of  rc and Bi leads to a decrease in the transient average temperature of all the layers.The proposed approach is verified by comparing the results with the finite difference solutions.The study shows that the proposed improved lumped model approach can be utilized as an effective analytical tool for the transient temperature behavior of multilayer spherical media.

Figure 1 :
Figure 1: Schematic illustration of cross-sectional view of -layer spherical media.

Figure 2 :
Figure 2: Comparison between the dimensionless transient temperature predicted by lumped model and the FD solution for each layer in the TRISO particle.

Figure 3 :
Figure 3: Temporal variations of the average temperatures predicted by the lumped model for different Φ 1 .

Figure 4 :
Figure 4: Temporal variations of the average temperatures predicted by the lumped model for different  rc .

Figure 5 :
Figure 5: Temporal variations of the average temperatures predicted by the lumped model for different Bi.