A Gaussian process regression accelerated multiscale model for conduction-radiation heat transfer in periodic composite materials with temperature-dependent thermal properties

Prediction of the coupled conduction-radiation heat transfer in composite materials with periodic structure is important in high-temperature applications of the materials. The temperature dependence of thermal properties complicates the problem. In this work, a multiscale model is proposed for the conduction-radiation heat transfer in periodic composite materials with temperature-dependent thermal properties. Homogenization analysis of the coupled conduction and radiative transfer equations is conducted, in which the temperature dependence of thermal properties is considered. Both the macroscopic homogenized equations and the local unit cell problems are derived. It is proved that the macroscopic average temperature can be used in the unit cell problems for the first-order corrections of the temperature and radiative intensity, and the calculations of effective thermal properties. The temperature dependence of thermal properties only influences the higher-order corrections. A multiscale numerical method is proposed based on the analysis. The Gaussian process (GP) regression is coupled into the multiscale algorithm to build a correlation between thermal properties and temperature for the macroscale iterations and prevent the repetitive solving of unit cell problems. The GP model is updated by additional solutions of unit cell problems during the iteration according to a variance threshold. Numerical simulations of conduction-radiation heat transfer in composite with isotropic and anisotropic periodic structures are used to validate the proposed multiscale model. It is found that the accuracy and efficiency of the multiscale method can be guaranteed by using a proper variance threshold for the GP model. The multiscale model can provide both the average temperature and radiative intensity fields and their detailed fluctuations due to the local structures.


Introduction
Prediction of the heat transfer in composite materials with periodic structures is important in many applications, from the insulation of buildings to the thermal protection of space vehicles [1][2][3]. Among a variety of prediction methods, homogenization techniques combine both macroscopic and microscopic solutions and can provide multiscale results in periodic composites [4,5]. The temperature field is divided into a macroscopic temperature field and local perturbations due to microscopic structures. The asymptotic analysis derives macroscopic governing equations for the whole material and governing equations in a representative unit cell. It also establishes the relation between macroscopic effective thermal properties and results of unit cell problems. The homogenization methods have been widely used in calculations of heat transfer in periodic composites [6][7][8], and are extended to the coupled conduction convection and radiation heat transfer problems [9][10][11][12][13][14][15][16]. The radiative transfer equation is also included in recent research [17].
In many works, the thermal properties of materials are assumed to be constant in the analysis. In practice, the thermal properties of real materials often depend on temperature. Therefore, the heat conduction in composite materials with temperature-dependent thermal conductivities is studied in some research. For example, Muliana and Kim [18] proposed a two-scale homogenization framework for laminated composites. The thermal conductivities of the components in the micromechanical model were temperature-dependent. Chung et al. [19] applied asymptotic expansion homogenization for composite materials with nonlinear properties. The thermal conductivities of the components in the unit cell problem were determined by the macroscopic average temperature or the non-uniform local temperature. Numerical tests demonstrated that the differences between the two approaches were insignificant. Zhai et al. [20] studied the transient heat conduction in composite materials with temperature-dependent thermal properties. The macroscopic average temperature was also used in the determination of the thermal properties in their unit cell problem. However, detailed asymptotic analysis of the heat conduction equation with temperature-dependent thermal properties is still needed.
Another problem associated with the multiscale model with temperature-dependent thermal properties is that the computational cost will be expensive. The effective thermal properties for the macroscopic model are provided by the solutions of the microscopic model. Because the temperature influences the thermal properties in the microscopic model, the microscopic model should be solved repetitively at each macroscopic iteration point, and the computational cost is unacceptable. This problem can be solved by reducing the order of the microscopic model [21]. For example, Monteiro et al. [22] used the proper orthogonal decomposition to reduce the order of the micro problem and simplify the computation. Another convenient approach is to build lists, fitted equations or surrogate models for the thermal properties to substitute the repetitive microscopic calculations. These surrogate models can be built before the macroscopic calculation based on preliminary microscopic calculations. They can be also built "on-the-fly", which means the surrogate models are updated during the calculations of macroscopic models [23,24]. In recent years, the rapid developing artificial intelligence and machine learning models provide tools for the multiscale methods. The machine learning models have been used in the hybrid atomistic-continuum simulations to substitute the microscopic molecular dynamics and complement the macroscopic continuum model [25,26]. Stephenson et al. [27] employed the Gaussian process (GP) regression model in their hybrid model. The advantage of the GP model is that it estimates the errors of the predictions simultaneously [28]. A threshold of the error can be used to automatically indicate when the update of the GP model and the additional microscopic calculations are needed. Therefore, it is promising to combine the machine learning model with the homogenization model to improve its efficiency for problems with temperature-dependent thermal properties.
In this work, a homogenization method will be proposed for the coupled heat conduction and radiative transfer equations with temperature-dependent properties, and the GP regression model will be employed to accelerate the multiscale simulations. The rest of the paper is organized as follows. In Section 2, an asymptotic analysis will be conducted on the heat conduction equation and radiative transfer equation. The temperature dependence of the thermal properties will be considered in the analysis and a homogenization model will be established. In Section 3, the GP regression model will be introduced first. An accelerated multiscale algorithm based on the GP regression model will be proposed. Some numerical examples are used to validate the proposed multiscale model in Section 4. Finally, conclusions will be presented in Section 5.

Homogenization of conduction and radiative transfer equations
The coupled heat conduction equation and radiative transfer equation are as follows [29]: is the coordinate and Ω = (Ω 1 , Ω 2 , Ω 3 ) is the direction of the radiative intensity. The T ε and I ε are the temperature and the radiative intensity. The heat conductivity k ε ij , the extinction coefficient β ε , the absorption coefficient α ε and the scattering coefficient σ ε are all functions of the coordinate and the temperature. The σ B is the Stefan-Boltzmann constant. The Einstein summation convention is used in this paper, and the repetitive indices i and j imply summations over all the components.
The homogenization of the Eqs. (1) and (2) is similar to the analysis in Ref. [17]. A local coordinate y = x/ε in a unit cell is introduced, as shown in Fig. 1. The ε can be regarded as the ratio between the characteristic length of the unit cell and the characteristic length of the whole material or computational domain. Therefore, ε is a small parameter when the material contains a large number of unit cells.
Then, the temperature and the radiative intensity are expanded as: where the T 1 and T 2 are periodic functions in a unit cell. It is also assumed that the volumetric averages of T 1 , T 2 and I 1 in a unit cell are zero. The spatial derivative becomes: The thermal properties k ε ij , α ε , β ε and σ ε are functions of both y and T ε . In the present work, the Taylor expansions are used and the thermal properties are expressed as: where φ represents any of the physical properties.
The Eqs. (3)-(6) are then submitted into Eqs. (1)-(2) to obtain the equations in different orders of ε. The detailed analysis can be found in Ref. [17], and the difference is that more terms are introduced into the equations by Eq. (6). However, it will be demonstrated that the temperature-dependent properties do not influence the macroscopic average equations and the first-order corrections.
The equation in the order of ε −2 still provides: The heat conduction equation in the order of ε −1 becomes: where the last term on the left-hand side is related to the variable thermal conductivity. Because ∂T 0 ∂y j = 0 , the first and last terms on the left-hand side can be eliminated. Therefore, the equation for unit cell problem is the same as that for problems with constant thermal properties, which is [17]: with the following form of T 1 : The periodic boundary condition is still used for N α and the volumetric average of N α in a unit cell should be zero. Eq. (9) demonstrates that the thermal conductivity in the unit cell problem is determined by the macroscopic average temperature T 0 . This analysis is also justified by the simulations in Ref. [19,20], where the T 0 is used in the unit cell problems because of its simplicity.
As for the radiative transfer equation, the ε −1 order equation still provides: In the order of ε 0 , the heat conduction equation gives: Here ∂T 0 ∂y j = 0 is already applied to eliminate a vanished term. The last term on the left-hand side is the additional term due to the thermal conductivity. The homogenized macroscopic equation is derived by integrating Eq. (12) in a unit cell. Because of the periodic properties of the variables in the unit cell, the terms started with ∂[] ∂y i will be eliminated.
Therefore, the temperature-dependent thermal properties still have no influence on the macroscopic heat conduction equation, and the homogenized equation is: The effective thermal conductivity K ij is given by: where |Y| is the volume of the unit cell, and the volumetric average value of a variable φ is defined as: The equations are the same as the equations derived with the constant thermal properties [17]. Therefore, the analysis further demonstrates that the macroscopic average temperature can be used in the unit cell problem to calculate the effective thermal properties.
However, the equations for the calculation of T 2 are more complex because of the additional terms and the dependence of T 0 on x. The spatial derivatives of k ε ij and N β are: Then, Eq. (12) can be reformed into: where all the properties are determined at T 0 . Equation (13) can be also rewritten as: The unit cell problems related to T 2 are derived by subtracting Eq. (18) from Eq. (17). According to the form of Eq. (17), the following ansatz for T 2 can be used: The corresponding governing equations for the M αβ , P αβ and C are: Compared with the analysis for problems with constant thermal properties in Ref. [17], it can be found that the equations for M αβ and C are the same, and the P αβ is the additional functions that should be solved. It is also easy to prove that the integrals of the right-hand sides of Eqs. (20)- (22) are all zero according to the periodic property of k ε and N α in a unit cell and the definition of K αβ in Eq. (14). The periodic boundary conditions can be still used for M αβ , P αβ and C in a unit cell. The previous research shows that the influence of T 2 on the accuracy of the temperature field is insignificant [17]. Because the present work mainly focuses on the computation of multiscale heat transfer problems with temperature-dependent thermal conductivities, the solution of T 2 is not considered and only T 0 and T 1 are calculated in this work for the temperature.
Finally, because the partial derivatives of thermal properties of the Eq. (6) are in the order of ε, they have no influence on the radiative transfer equations in the order of ε −1 and ε 0 . Therefore, the equations for I 0 and I 1 are the same as those in Ref. [17]: It should be mentioned that a variety of boundary conditions of I 1 can be used for Eq. (24) in a unit cell. In the previous work, it is found that the periodic boundary condition will overpredict the fluctuations in the unit cell and lead to large errors [17]. It is suggested that the I 1 is solved in the whole region with the following boundary condition: where n is the outward unit normal and Γ is the boundary of the whole domain. The I 1 solved by Eq. (25) is denoted by I 1,all in this work. However, this treatment of I 1 still has two drawbacks. Firstly, the error of the radiation tends to accumulate from the incoming boundary to the outgoing boundary, and the fluctuations will be amplified. Secondly, it is inconvenient to obtain a local radiation fluctuation by solving I 1 in the entire region. In this work, the following I 1 solved in an individual unit cell is used: where Γ Y is the boundary of a unit cell. The I 1 solved by Eq. (26) is denoted by I 1,cell in this work. The boundary condition in Eq. (26) means that the incoming I 1 is zero in each unit cell, which is different from the periodic boundary condition. It is found that I 1,cell tends to underpredict the fluctuations. Therefore, the average of I 1,all and I 1,cell is also used as a boundary condition, which is calculated by: In summary, the influence of the temperature-dependent thermal properties only appears in the equations for T 2 . The macroscopic homogenized equations for T 0 and I 0 and the unit cell problems for T 1 and I 1 are the same as the equations for problems with constant thermal properties. The macroscopic average temperature T 0 is used in these equations to determine the thermal properties. All the governing equations are summarized in Table 1 for convenience.

Gaussian process regression accelerated multiscale algorithm
According to the analysis in Section 2, a multiscale computational procedure can be proposed for the conduction-radiation heat transfer in periodic composite materials. A mesh for the macroscopic simulations and a mesh for the unit cell problems are firstly established. Then, the Eqs. (13) and (23) are solved to obtain the macroscopic T 0 and I 0 . The effective thermal conductivities and radiative transfer coefficients are calculated from the solutions of the unit cell problems. Because the properties in the unit cell problems depend on the temperature, the above calculations of the effective properties should be employed in each iteration on all the mesh nodes. This procedure is unacceptable because the repetitive calculations of the unit cell problems under different temperatures will consume large computational resources.
Therefore, in the present work, the Gaussian process regression will be used as the surrogate model to obtain the correlation between the effective thermal conductivity and the temperature. In the rest of this section, the GP regression will be briefly introduced in Section 3.1, and the multiscale computational procedure will be described in Section 3.2.

Gaussian process regression
Details of the GP regression model can be found in Ref. [27,28], and only the essential algorithms are described in this section. It is assumed that n thermal conductivity tensors K ij,k , k = 1 ~ n, have been calculated from the unit cell problem before the regression, which are corresponding to the temperatures T k . The vectors K ij and T are defined as Table 1 The summary of the equations for the homogenization of conduction and radiative transfer equations 1 , K ij,2 , …, K ij,n ) T and T = (T 1 , T 2 , …, T n ) T . In the GP model, the prior mean function is chosen as zero. The variance of the noise of K ij is denoted by σ 2 n . The σ n is chosen as the residual of the N α , which is 10 -7 in this work. The squared exponential covariance function is used in this paper:

Variables Governing equations
where σ f and l are two hyperparameters. Then, the probability distribution of K ij * at temperature T * is a Gaussian distribution with mean value [28]: and the variance: where C * is an n × 1 vector with entries C * k = C(T * , T k ) and C is an n × n matrix with entries C km = C(T k , T m ).
During the simulation, when the C and K ij are calculated based on the existing data, the K ij * at a new temperature T * can be predicted by Eq. (29). A variance σ 2 * is also given by Eq. (30). If the σ * exceeds a threshold σ t which is specified before the simulation, a new unit cell problem of N α will be called to calculate the K ij at the new temperature.
The new data will be added to the existing data set and the C and K ij are updated. In addition, the σ * is not relevant to K ij , as demonstrated by Eq. (30). Therefore, the threshold is the same for all the components of thermal conductivity tensor.
In practice, the two hyperparameters in the GP model need to be determined. A possible method is to use the maximum likelihood estimation (MLE) as mentioned in Ref. [27]. In the present work, we simply choose l = T h -T c and σ f = 1. The reason can be explained by the following test. The data is the K 12 calculated from a unit cell with elliptical particle with temperature changing from 400 to 1400 K. The K 12 for 400 K, 800 K, 1200 K and 1400 K are used to train a GP model. The MLE gives an optimized set of (σ f , l) as (0.0075, 370.6). However, the K 12 predicted with the optimized hyperparameters have large deviations compared with the real data, as shown in Fig. 2. On the contrary, the GP model with (σ f , l) = (1, 1000) coincides better with the calculated data. The reason is that although the optimized hyperparameters by MLE are more consistent with the training data, the training data could not reflect the trend of the real data. It should be mentioned that although different hyperparameters provide different predictions, the parameter that controls the accuracy of the multiscale simulation is the threshold σ t . In the present work, the (σ f , l) = (1, 400) is used and the influence of σ t will be discussed.

Multiscale computational procedure
The computational procedure for multiscale simulation with GP regression is given in Fig. 3. The procedure contains the following steps.
(1) Build meshes and initialize the macroscopic problems and the unit cell problems.
Select a set of initial temperatures. Solve Eq. (9) for N α in the unit cell and calculate the effective thermal conductivity and radiative coefficients by Eqs. (14) and (15).
(3) Build an initial GP model for the thermal properties by calculating the matrix C(T, T) + σ 2 n I −1 in Eqs. (29) and (30).
(4) The macroscopic homogenized Eqs. (13) and (23) are iterated, while the thermal properties are calculated from the GP model. Equation (30) is firstly used to obtain the σ * according to the current temperature of the grid node. If σ * is smaller than the threshold σ t , Eq. (29) is used to calculate the thermal properties. Otherwise, the unit cell problem is solved to obtain the new thermal properties under the temperature. The GP model is updated based on the enriched database of the thermal properties. (5) Equations (13) and (23) are iterated until convergent results of T 0 and I 0 are obtained. (6) Finally, the T 1 and I 1 in each unit cell can be obtained by Eqs. (9), (10) and (24). The first-order approximations for T ε and I ε can be reconstructed by: In this work, the Cartesian meshes are used. The finite volume method and discrete ordinate method are employed to solve the heat conduction and radiative transfer equations, respectively. Details of the numerical methods can be found in Refs. [17,30]. In addition, it should be mentioned that the value of the small parameter ε used in the computations has no influence on the final results. The reason is that the solutions of the equations in the unit cell, such as Eq. (9) and Eq. (24), are inversely proportional to the ε in the coordinate y. Thus, the ε is cancelled out in Eqs. (31) and (32) when calculating the final results. More detailed explanations can be found in Section 3 of Ref. [17]. In this (31)  The range of the temperature T in Eq. (33) for the correlations in Eqs. (34) and (35) is 300 K ~ 1400 K.
The first example is a 1 mm × 0.1 mm region with ten 0.1 mm × 0.1 mm unit cells [17], as shown in Fig. 4. The TiO 2 particle is assumed to be a circle with a diameter of 20 μm. The temperature T h on the left boundary is 900 K, 1100 K or 1300 K. The temperature T c on the right boundary is specified as 400 K. The boundaries are assumed to be black boundaries for the radiation. Periodic boundary conditions are used on the other boundaries. For the fully-resolved simulations, the size of the mesh for a unit cell is 400 × 400 and that for the whole domain is 4000 × 400. The grid independence for effective thermal conductivities can be found in Ref. [17]. As for the homogenized equations, the grid size is 1000 × 100.
The criterion for convergence is that the relative error between two iterations is less than 10 -8 . Because the grid size is relatively large and the thermal conductivities are temperature-dependent, the convergence for fully-resolved simulation is slow. Therefore, in this work, the temperature results of the macroscopic homogenized equations are used as the initial temperatures for the fully-resolved simulations to speed up the convergence. The relative errors of the temperature and radiative intensity fields of the multiscale simulations are defined as: where the subscripts i and j are the indices for 2D grid nodes and k is the index for discrete directions of radiative intensities. The superscripts M and F represent the results of the multiscale simulations and the fully-resolved simulations. The reference temperature is the average of the two boundary temperatures: T ref = (T h + T c ) 2 , and the reference radiative intensity is Taking the problem with 1300 K left boundary temperature as an example, the unit cell problem is firstly solved under 6 temperatures: 400 K, 600 K, 800 K, 1000 K, 1200 K and 1300 K. An initial GP model for the effective thermal conductivities is established according to the data on the 6 temperatures. When the σ * of the GP model on a new temperature is larger than the threshold σ t , the unit cell problem is solved, the data set is updated and a new GP model is built. The influence of the value of σ t is considered by setting σ t as 10 -2 to 10 -6 . The results for the multiscale simulations are given in Table 2. Evidently, when σ t decreases, the number of temperature points on which the effective thermal conductivity is calculated increases. The CPU time increases because more unit Fig. 4 Computational domain for Example 1 cell calculations are conducted and the dimension of the matrices in Eqs. (29) and (30) increases. As a comparison, the CPU time for the fully-resolved simulation is 69.2 h, although the solution of the homogenized equation is used as the initial temperature. Thus, the multiscale method increases the efficiency of the simulations. Meanwhile, it can be found that σ t = 10 -4 is sufficient for the simulation because a smaller threshold cannot further increase the accuracy. Therefore, a proper threshold can be chosen to guarantee both the accuracy and efficiency of the multiscale model, and σ t = 10 -4 is used in the rest of the work. The temperature and incident radiation fields of the fully-resolved simulation and the multiscale simulation with σ t = 10 -4 are shown in Fig. 5. The incident radiation in Fig. 5 is defined as G = k w k I k , where w k and I k are the weights and discrete radiative intensities in the discrete ordinate method. Both the I 1,cell in Eq. (26) and the I 1,ave in Eq. (27)

Fig. 5
The temperature and incident radiation fields of the fully-resolved simulation and the multiscale simulation are used to reconstruct the radiative intensity fields. It can be seen that the temperature and incident radiation fields of both the simulations coincide well with each other. The multiscale model can provide detailed small scale fluctuations of temperature and radiative intensity fields. The relative errors of the incident radiation are 12.47% for I 0 , 10.72% for I ε with I 1,cell , and 5.90% for I ε with I 1,ave . Fig. 5 (b) also demonstrates that the I 1,cell underpredicts the fluctuations and I 1,ave provides a better result. The 900 K and 1100 K are also used as the left temperature boundary conditions to show the influence of macroscopic temperature gradient on the multiscale simulations. The errors of temperature and radiative intensity are given in Table 3. The results demonstrate that the multiscale model can provide accurate results of the coupled conduction-radiation heat transfer in composites. The temperature and incident radiation profiles along horizontal line y = 0.05 mm are also shown in Fig. 6. It can be seen that the fluctuations of temperature and radiation induced by the particle are reproduced by the multiscale model.
Finally, the influence of radiation on the accuracy of the multiscale method is considered. The 1300 K and 400 K are used as the temperature boundary conditions. The grid size is still 4000 × 400. The absorption and scattering coefficients are multiplied by a factor to tune the effects of radiation. The relative errors of temperature and radiative intensity fields for different multipliers are given in Table 4. Although the influence on the error of temperature is insignificant, the error of radiation intensity increases with the increase of absorption and scattering coefficients. The reason can be found in Fig. 7 and Fig. 8, where the profiles of incident radiations with different multipliers are presented. The α 2 = 72080 m −1 corresponds to the multiplier 10 in Fig. 8. The profiles have similar      Table 5 and Fig. 8. It can be seen that the errors of the radiation decrease as the values of the absorption coefficients of the two components become closer and the fluctuations become weaker.

Example 2
In this section, the coupled conduction and radiation heat transfer with anisotropic structures are simulated [17]. The structures are shown in Fig. 9, which are composed of 7 × 7 unit cells. The size of the unit cell is 0.1 mm and the size of the domain is 0.7 mm. The lengths of the major and minor axes of the elliptical particles are 94.9 μm and 9.49 μm. The particles are horizontally aligned or 45° tilted. The properties of the materials are the same as in Example 1. The temperature on the left boundary is specified as T h = 1300 K and the temperature on the other boundaries is T c = 400 K. The grid size of a unit cell is 400 × 400, and that of the fully-resolved simulation is 2800 × 2800. The grid size of the macroscopic homogenized simulation is 700 × 700. The temperature and incident radiation fields for both the fully-resolved simulations and the multiscale simulations are given in Fig. 10 and Fig. 11. It can be seen that the multiscale results can reproduce the local temperature and radiation fluctuations due to the elliptical particles (Table 6). In addition, the CPU times of the fully-resolved simulations are 596.1 h and 500.2 h for the horizontally aligned and tilted particles, respectively. The CPU times of the multiscale simulations are 11.7 h and 13.7 h.

Conclusions
In the above research, the homogenization of the coupled conduction and radiative transfer equations with temperature-dependent thermal properties has been conducted. The Taylor expansions of the thermal properties about the average temperature are used in the analysis. Both the homogenized equations and the unit cell problems are derived. It is proved that the macroscopic average temperature can be used in the unit cell problems to determine the effective thermal properties. The homogenized equations for T 0 , I 0 and the equations for T 1 , I 1 are the same as the equations in the analysis with constant thermal properties. The effect of the temperature-dependent thermal properties only occurs in the higher-order corrections such as the equation for T 2 . Based on the homogenization analysis, the multiscale numerical method has been proposed. The cell problem is solved to provide effective thermal properties for the macroscale iterations. In order to prevent the repetitive calculation of the cell problem, the GP regression is introduced to build a correlation between the effective thermal properties and the temperature. The GP model is updated during the iteration according to the variance of the model.
By comparing with the results of fully-resolved simulations of conduction-radiation heat transfer in composite with isotropic and anisotropic periodic structures, the proposed multiscale method is validated. By using a proper variance threshold for the GP Fig. 10 Temperature and incident radiation fields for horizontally aligned particles model, the accuracy and efficiency of the multiscale method can be guaranteed. The multiscale method can provide both the average temperature and radiative intensity fields and their fluctuations due to the local structures. It is found that the I 1,ave can provide better correction of the radiative intensity, which is the average of the I 1,all solved in the whole domain and the I 1,cell solved in the individual cell with zero boundary incoming. It is also found that the error of the radiation increases with the magnitude of the absorption and scattering coefficients and the difference between the coefficients of different components. The proposed multiscale model can be used in further studies on the high-temperature heat transfer in composite materials with periodic structures. It should be mentioned that the analysis in this work is for the interior of composites. The  boundary layer problems should be studied in the future for the correction of boundary conditions. More detailed analysis of the homogenization of the radiative transfer equation is also needed to obtain a clearer boundary condition for I 1 [31,32].