Reproducing Kernel Particle Method for Radiative Heat Transfer in 1 D Participating Media

The reproducing kernel particle method (RKPM), which is a Lagrangian meshless method, is employed for the calculation of radiative heat transfer in participating media. In the present method, for each discrete particle (i.e., spatial node) within a local support domain, the approximate formulas of the radiative intensity and its derivatives are constructed by the reproducing kernel interpolation function, and the residual function is obtained when these parameters are substituted into the radiative transfer equation. Then the least-squares point collocation technique (LSPCT) is introduced by minimizing the summation of residual function. Five test cases are considered and quantified to verify the meshless method, including isotropic scattering medium, first-order forward scattering medium, pure absorbing medium, absorbing scattering medium, and absorbing, scattering emitting medium. The results are in good agreement with the benchmark methods, showing the reproducing kernel particle method is an efficient, accurate, and stable method for the calculation of radiative transfer in participating media.


Introduction
Radiative heat transfer among absorbing, scattering, and emitting medium plays an important role in many engineering application, that is, glass fabrication, laser pulse heating, and diesel fuel droplets burners.As a useful tool, numerical simulation is adopted widely in the study of radiative heat transfer, in which some methods have been developed to solve the radiative heat transfer in participating medium, such as the Monte Carlo method, the zonal method, the discrete-ordinates method (DOM), the finite volume method (FVM), and the finite element method (FEM).However, these traditional methods severely depend on the predefined mesh quality for complex geometry shape and remeshing is needed while large deformation occurs.As a result, some meshless methods are proposed to avoid the drawback.
Sadat [1] presented a moving least-squares collocation method to solve the radiative transfer equation, and the results showed that, for small absorption coefficient, the primitive variables formulation is unstable while the even parity formulation is always stable and accurate.Liu employed the meshless local Petrov-Galerkin (MLPG) method [2][3][4] and the least-squares collocation meshless (LSCM) method [5][6][7] to solve radiative heat transfer in gradient refractive index medium and coupled heat transfer in a participating medium.Wang et al. [8] presented a moving least squares approximation meshless method to solve the radiative transfer equation in complex 2D and 3D geometries.Zhang et al. [9,10] proposed the Nature Element Method (NEM) to simulate the radiative heat transfer in a two-dimensional enclosure with absorbing, emitting, and isotropically scattering medium inside.
The reproducing kernel particle method (RKPM) [11] was developed from the smoothed particle hydrodynamics (SPH) method [12].Both of the methods are Lagrangian meshless methods.Compared with other meshless methods, where the spatial node is only interpolated node, the Lagrangian meshless method has itself advantage that the discrete particle carries volume, density, and mass of material and can move under effect of external force and inner interaction.The discrete particle represents both the approximate node and material characteristics at same time, and it makes the Lagrangian meshless method more attractive.However, SPH has a drawback, where the shape function does not satisfy normalization when the discrete particle distribution is nonuniform, and this also happens on the boundary when the discrete particle distribution is uniform.In RKPM, the shape function was improved by introducing a continuous reproducing kernel.So the RKPM can work well over the whole domain.It has wide application in many engineering fields, such as, rolling plane strain problem [13], bucking analysis of thin plates [14], large deformation nonlinear elastic problems [15][16][17][18], metal forming problem [19][20][21][22], elastic-plastic problems [23,24], convection-diffusion problem [25][26][27], heat conduction problems [28][29][30], and fragment-impact problem [31,32].But up to now, to the best of our knowledge, there is still lack of literature on the application of RKPM for radiative heat transfer.
The objectives of this paper are to introduce the RKPM for solving the problem of radiative heat transfer in participating medium, in which the RKPM is employed to construct the interpolation function of the radiative intensity and its derivative and the least-squares point collocation technique (LSPCT) is employed to obtain the numerical solution by minimizing the summation of residuals of all discrete particles and determine the optimized scaling factor value and number of discrete particles in calculated process.In Section 2, the principle of RKPM is described briefly and the discretization of radiative transfer equation is constructed.In Section 3, five numerical examples are introduced to verify the performance of the method.Finally, conclusions and future work on the RKPM are given in Section 4.

Reproducing Kernel Particle Method (RKPM).
The RKPM has a similar approximation formula with the SPH method; therefore the integral physical field function (  ) can be written as follows [33]: where ,   are the spatial coordinates; ℎ is the radius of the support domains of the corrected kernel function, and ℎ = Δ, in which  is a scaling factor and takes an empirical value, where Δ is the distance of two adjacent discrete particles in the -direction; () and (  ) are the physical field functions at spatial coordinates  and   , respectively; Ω is the supported domain of spatial coordinate ; (ℎ,  −   ) is the corrected kernel function and the product of correction function (,   ) and kernel function (ℎ,  −   ).
For the RKPM, the interpolation function is constructed on the discrete particle (i.e., spatial node, which contains partial volume of studied domain) in the whole domain, and the discretization of (1) could be written as follows: where   and   are spatial location of discrete particle and discrete particle in the support domain of the th discrete particle, respectively;   is the total number of the discrete particle ; Δ  is the length for a one-dimensional problem, or the area for a two-dimensional problem, or the volume for a three-dimensional problem;   is the shape function of the discrete particle.
The -derivative of the approximation formula can be written as where  , is the derivative of the shape function for th discrete particle.

Discrete Ordinate Equation.
Discrete ordinate formula of one-dimensional radiative transfer equation is as follows [35]: where  and  are one discrete direction and total discrete direction of radiation, respectively;   is the radiative intensity of the th direction in medium;   is the black body radiative intensity in medium;   and   are absorbing coefficient and scattering coefficient in medium, respectively;   is the cosine value of th direction on -direction;   is the solid angle weight of th direction.Φ    is the scattering phase function;  is the refractive index of medium.
The radiative boundaries are considered to be opaque and diffuse as well as gray.The radiative intensity of the wall can be written as follows: where  represents the wall;   is the emissivity of the wall.n  is the outside unit normal vector of the wall.s   is the unit vector on   direction.
For ( 5), the forward scattering term is moved from the right side to the left side.Equation ( 5) can be rewritten as

Discretization and Numerical Implementation.
The radiative intensity and its derivatives are approximated by RKPM and substituted into the radiative transfer equation to obtain the residual function of the radiative transfer problem.The least squares point collocation technique (LSPCT) is adopted to solve the equations.The basic idea of LSPCT is to obtain the numerical simulating solution by minimizing the summation of residual functions in the whole domain.The numerical implementation is described in the following section.
According to (4a) and (4b), the approximations of the radiative intensity and its derivatives were constructed on the th discrete direction at the th discrete particle, shown as follows: where    is the radiative intensity of the th discrete particle on the th discrete direction.
The residual function can be obtained by substituting (8a) and (8b) into (7), thus where (  ) is the residual function of physical field function for each discrete particle . are the total numbers of discrete particles in the whole domain: Then, the summation of the residual function for all the discrete particles can be obtained: Minimizing the function ((  )) with respect to the radiative intensity   (  ), the derivative can be written as Moreover, the linear matrix equation can be expressed as where For each discrete particle  on the inflow boundary, the radiative intensity is given by ( 6), and the boundary condition can be directly imposed as follows: where   is the Dirac delta function.
The calculation steps for radiative transfer are as follows: firstly, the shape function of RKPM is calculated according to the current distribution of discrete particle and substituted into (10a) and (10b); secondly, the matrix coefficient of ( 13) is obtained by substituting (10a) and (10b) into (14a) and (14b); finally, the radiative intensity for all the discrete particles is calculated and used to solve (13).Due to the scattering effect, the radiative intensity would reiterate the process for all directions until the convergence condition is satisfied.(b) FVM of [37] Figure 1: Effects of discrete particle number and optical thickness for relative error between the numerical results from this meshless method and other different numerical methods in Case 1.

Results and Discussion
3.1.Radiative Transfer in Pure Scattering Medium.In this section, a one-dimensional radiative transfer is considered in pure scattering gray medium, for which the surfaces are black body, and the temperature of the left and the right surface is 1000 K (  1 ), 0 K (  2 ), respectively.The incident radiative heat flux   of the right surface is defined as (16), and the dimensionless reflection radiative heat flux is defined as q   = 1 −   /( 4  1 ), which is analysed in the two different scattering mediums.The scattering phase functions are isotropic scattering phase function and first-order forward scattering phase function, respectively.The relative error is defined with (17) for the dimensionless reflection radiative heat flux to compare the results from this method and the reference: Case 1 (isotropic scattering medium).The scattering phase function expression is The relative error of the dimensionless radiative heat flux is defined as shown in (17), and   RKPM and   Ref. separately represent the dimensionless radiative heat flux is in present work and the corresponding data in [36,37], and the relative error is affected by the discrete particles number and optical thickness as shown in Figures 1(a) and 1(b).Through analyzing these two figures, we can know that the variational trend of relative error is irregular when the discrete particles number is 11.When the discrete particle number is larger than 11, the relative error is decreasing with the increase of discrete particles number.And the relative error remains at the same value when the discrete particle number exceeds 81.Therefore, the optimal number of discrete particles is recommended 81 in the following calculation.
Here, we can also find that the relative error is decreasing with the increase of optical thickness as shown in Figures 1(a) and 1(b).And the value of relative error is very small and its maximum is only 0.537%.So this meshless method has a high accuracy when computing the radiative transfer for the optical thickness ranging from 1 to 10.
The dimensionless radiative heat flux varies with the optical thickness and scaling factor is shown in Table 1.Good agreement is observed for the results from the present method, the theory solution from [36] and the numerical solution from [37].Since the impact on the flux is small while changing the scaling factor from 1.01 to 1.10, the scaling factor is recommended within this range to calculate the radiative transfer in 1D isotropic scattering medium.
The time required for computation in Case 1 is shown in Table 2.The computation is implemented on a personal computer with an Intel Celeron Dual-Core 1.8 GHz Processor.Increasing scaling factor value does not increase the time required for computation much more, and increasing the number of discrete particle and optical thickness increases the time required for computation much more.
Case 2 (first-order forward scattering medium).For the first-order forward scattering medium, the scattering phase function can be written as Φ    = 1 + cos   cos    .
Figures 2(a), 2(b), and 2(c) show the relative error of the dimensionless radiative heat flux, and the values in this function come from the present work and the corresponding data in [38][39][40].As the same to Case 1, from the three figures, (c) FVM of [40] Figure 2: Effects of discrete particle number and optical thickness for relative error between the numerical results from this meshless method and other different numerical methods in Case 2. the numerical results are not good when the discrete particles number is 11.The relative error is decreasing with the increase of discrete particles number, when the discrete particle number is larger than 11.Then the error remains at the constant value when the discrete particle number exceeds 81.However, relative error has different variational trends in the three figures.In Figure 2(a), the variational trend increases at first and then decreases with the increase of optical thickness.In Figures 2(b) and 2(c), the variational trend decreases with the increase of optical thickness.The reason is that the results are obtained from the numerical method and depended on the feature of the algorithm.And the maximum value of relative error is small (less than 0.923%).It shows that this meshless method has good value to be applied in radiative transfer region.Table 3 shows the effects of the scaling factor and the optical thickness.The results using the present method are in good agreement with the data of [38][39][40].Comparing the results in the isotropic scattering medium and in the first-order forward scattering medium, it is found that their variation trends are very similar.That means there should be the same principle for the two different kinds of scattering medium in terms of selecting the number of discrete particle and the value of the scaling factor.
The time required for computation in Case 2 is shown in Table 4.The computation is implemented on the same computer as Case 1.And the time required for computation has the same trend with Case 1.
The results obtained from Cases 1 and 2 show that the RKPM can simulate radiative transfer in pure scattering medium, which contains isotropic and anisotropic medium.In the two cases, when the optical thickness of media ranges from 1 to 10, the simulative results are accurate comparing with the data in the corresponding reference.This verifies that the RKPM is applicative and accurate for solving radiative transfer in pure scattering medium.

Radiative Heat Transfer in Absorbing Medium
Case 3 (radiative transfer in pure absorbing medium).The radiative heat transfer in a one-dimensional pure absorbing medium is considered in this section.Both boundary walls are regarded as the black body, and the temperature of side walls is  1 (left side) and  2 (right side), respectively.The dimensionless temperature and dimensionless radiative heat flux density are defined as follows: In this study, the optical thickness value is set as 100 to replace the infinity.The discrete particles number, scaling factor of smooth kernel function, and iterative convergence error are 81, 1.01, and 1 × 10 −6 , respectively.Figure 3 shows that the dimensionless temperature distribution in the absorbing medium with different optical thicknesses.By comparison, a good agreement was observed between the data obtained by present method and the corresponding data from [41].Table 5 and Figure 4 show the effect of optical thickness on the dimensionless radiative heat flux and the relative error, in which it is clear that the relative error decreases with the increasing optical thickness.The maximum relative error is 1.4% when the optical thickness is 0.1.
The time required for computation in Case 3 is shown in Table 6.The computation is implemented on the same computer as Case 1. Increasing optical thickness increases the time required for computation much more, and the time required is 277.812s when the optical thickness is 100.0.In Case 3, the RKPM solves radiative transfer in pure absorbing medium under different optical thickness.The simulative results are accurate comparing with the data in the reference [41].This verifies that the RKPM is applicative and accurate for solving radiative transfer in pure absorbing medium.

Radiative Transfer in Absorbing, Scattering, and Emitting
Medium.The results obtained from Cases 1, 2, and 3 show that 81 is taken as the number of discrete particles and 1.01 as the scaling factor value for the calculation of radiative transfer, which are used for Cases 4 and 5.
Case 4 (radiative transfer in absorbing and isotropic scattering medium).In this case, the optical thickness of absorbing and isotropic medium is 1, the Albedo's value is 0.5, and the initial temperature is 0 K.Both sides of the wall are black body, and the temperature of left surface is  1 , the right surface  2 .Both dimensionless incident radiation intensity and dimensionless radiative heat flux density are calculated in the medium.The incident radiation intensity and radiative heat flux density are written as follows: Case 5 (radiative heat transfer in absorbing scattering and emitting medium).In this case, the value of Albedo is 0.5, and the optical thickness is selected as 0.1, 1, and 10, respectively.Both sides of the walls are the black body, and the temperature of left surface is  1 , the right surface  2 .Distribution of the dimensionless incident radiative intensity, dimensionless radiation heat flux density, and dimensionless temperature is calculated.The temperature is defined as follows: It is obvious from Figures 5, 6, and 7 that the dimensionless incident radiative intensity, the dimensionless heat flux, and the dimensionless temperature by this method agree well with the corresponding data in [42] obtained by DOM, which means that the meshless method has good accuracy in solving the radiative transfer problems in absorbing isotropic scattering medium and absorbing isotropic scattering emitting medium.For Case 5, the distributions of dimensionless temperature change are calculated at three optical thicknesses,   = 0.1,   = 1, and   = 10 with this method, as shown in Figure 7.The gradient of dimensionless temperature curves    increases as the optical thicknesses become larger, which has same increasing trend for the data from RKPM and DOM.It means the radiation energy, coming from high temperature boundary, is not allowed to penetrate deeply into the medium when the optical thickness is large enough.And for a certain optical thickness, the curve in low temperature zone changes more steeply than it does in high temperature zone, especially for the large optical thickness.
In Cases 4 and 5, the RKPM simulates the radiative transfer in absorbing, scattering, and emitting medium under different optical thickness.The results are accurate comparing with the data in the corresponding reference.This verifies that the RKPM is applicative and accurate for solving radiative transfer in absorbing, scattering, and emitting medium.
The time required for computation in Cases 4 and 5 is shown in Table 7 by RKPM and DOM.The computation is implemented on the same computer as Case 1. Generally, the computation process of meshless method needs more time.The reason is that the meshless method spends more time to search adjacent particle and calculate the coefficient of polynomial.This is shown in Table 7.
From the above five cases, the simulative results are stable and accurate, which verify the application of RKPM for solving radiative transfer in participating medium under different optical thickness.

Conclusion
A reproducing kernel particle method based on the discreteordinate equation is extended to solve the radiative transfer problem in participating medium.The reproducing kernel interpolation function is used to construct the approximation formula of radiation and its derivative.For each discrete particle, the residual function is obtained by substituting approximation formula into radiative transfer equation.The least-squares point collocation technique is adopted to obtain the solution of the problem by minimizing the summation of residuals of all discrete particles.In the present work, the RKPM is a successful implementation to simulate radiative transfer in pure scattering and absorbing medium and in absorbing, scattering, and emitting medium.The results show that the RKPM is efficient, accurate, and stable and can be used for solving radiative heat transfer in participating medium.And the optimized scaling factor value and number of discrete particle are obtained in calculated process, which are between 1.01 and 1.10 and no less than 81, respectively.
Therefore, it can be achieved to simulate high temperature phase change of semitransparent material with RKPM.Based on the present work, the application of RKPM will be explored further to the radiative heat transfer in 2D and 3D semitransparent media and the combining conductionradiation heat transfer in high temperature semitransparent media with phase change.

Figure 4 :Figure 5 :
Figure4: The relative error of dimensionless radiative heat flux density between the simulated results from this meshless method and the data from[41] under different optical thickness.

Table 1 :
Effects of scaling factor value for dimensionless heat flow under discrete particle number 81 in Case 1.

Table 2 :
Time required for computation in Case 1.

Table 3 :
Effects of scaling factor value for dimensionless heat flow under discrete particle number 81 in Case 2.

Table 4 :
Time required for computation in Case 2.

Table 5 :
Dimensionless radiative heat flux density for different optical thickness.

Table 6 :
Time required for computation in Case 3 of  = 1.01 and  = 81.

Table 7 :
Time required for computation in Cases 4 and 5.