Analysis of Heat Transfer in the Case of Non-linear Hyperbolic Conduction Equation Coupled with Radiation or with Thermoelectricity

Simultaneously, transient conduction and/or radiation in participating media appears in many engineering systems and emerging technologies such as nanostructures, biological tissues, insulated foams, polymers, furnaces, heat pipes, combusting chambers and rocket propulsion, renewable energy applications using thermoelectric materials [1-3]. The amount and rate of heat, coupled to system properties and the surrounding govern the thermal response at both the microscopic and macroscopic scales. For transient heat flow during an extremely short period, the conduction process is well described by the hyperbolic heat conduction theory based on non-Fourier constitutive heat flux equations such as the Cattaneo – Vernotte heat flux equation, which implies a finite speed for the propagation of thermal perturbations [4-5]. In engineering practice, it is convenient to simulate the non-linear heat conduction behavior with commercial software such as ANSYS, FLUENT and TRNSYS [6-8]. However, it is interesting and suitable to develop semi analytical models for instance to make a thermodynamic analysis of a thermoelectric device [9], to estimate the maximum rate of heat pumping, to study the influence of one parameter, to use it for optimization or as an inverse model. Moreover, several modeling of the transient state is based on the convenient electrical analogy, Laplace transform and separation of the variables [10] finite difference/volume and the lattice Boltzmann (LBM) methods [11] to study thermoelectric self-cooling of devices. They are also convenient to solve the partial differential equation for a small time leg or for micro thermoelectric coolers [12].

At macroscopic scales, the time and spectral changes of radiative information are often considered as non-significant due to its very large propagation speed and gray medium assumption [13]. Hence, the radiation assessment involves an integro-differential radiative transfer equation (RTE) including three spatial and two angular variables. Furthermore, for present curvilinear shape, the differential form of the RTE coupled to the dependence of spatial/angular variables and the optically inhomogeneous complex media make the radiative transfer equation difficult to solve analytically except for some limiting cases [1, 2, 13]. Generally, the solutions of RTE found in the literatures are typically determined by numerical means and therefore, contain errors such as angular and spatial errors due to solution procedures. To mitigate the effects of error due to angular discretization, different approaches have been described in the literatures, including phase function renormalization, high order angular discretization, angular grid refinement [2, 11, 14-15]. However, the error due to spatial discretization depends to the considered approach around each grid node and to the selected method. In order to improve the accuracy of the solution, semi-analytical approach based on discrete ordinates has been developed recently by Ladeia et al. [16] for cylindrical media but limited to homogeneous and steady state transfer.

The present research is concerned: (1) firstly with the development of an analytical solution for radiation analysis in optically complex media of inhomogeneous graded index. The spherical harmonics method is then used for the development of an analytical radiative transfer solution. The non-linear hyperbolic heat conduction with temperature dependent thermal conductivity has been investigated through the lattice Boltzmann and finite volume methods; (2) and secondly with the complete one-dimensional transient hyperbolic heat transfer equation solved semi-analytically in order to obtain explicit expressions of the temperature and also of the heat flux in a thermoelectric layer. A quadrupole formulation of the hyperbolic partial differential equation taking the Thomson effect into account is then proposed for thin thermoelectric layer or medium considered at very short times.

The first solution approach consists to solve the radiation transfer using the spherical harmonics in order to build the solution of Eq. (1) from Lattice Boltzmann method. Secondly, a specific quadrupole method is developed to solve the energy equation relative to a thin layer applied to thermoelectricity.

3.1 A layered radiative solution

The inhomogeneity and angular redistribution of radiative intensity increase the difficulties to build the analytical solution of the radiation equation. In order to construct a semi-analytical solution, the media is assumed to be a set of sub-layers with constant factor 1/τ. In such case, the radiative intensity in a single layer can be written as

$μ (∂I_k)/∂τ+(1-μ^2 )  a/τ  (∂I_k)/∂μ+I_k-[1-ω] I_b (T)=ω/2 ∫_{-1}^1 Φ(μ'→μ) I_k (τ,μ' )dμ' $  (5)

where $ I_b=n^2 σ_B T^4/π$ is the blackbody intensity; $τ=κ_e r$ being the optical depth; $μ=cos⁡(θ)$ is the polar direction cosine; and $Φ$ is the scattering phase function. For the constant a=0 and 1, the equation pertains to that for the slab and the spherical enclosure, respectively.

The radiative problem in each layer is solved using the spherical harmonics or P_N method. Therefore, the radiative intensity is expanded as [2, 18]

$I_k (τ,μ)≈∑_{l=0}^NI_(l,k) (τ)P_l (μ)$    (6)

where $I_{l,k}$ are the radiative intensity moments to be determined and $P_l$ are the Legendre polynomials. Hence, applying recursive expression of $μP_l (μ)$ and combined to orthogonality condition of Legendre polynomials, the RTE is reduced to the matrix- vector form as

$A (dF^k)/dτ+(C+a/τ_k  B) F^k=D^k;F^k=[I_{0,k},I_{1,k}  ,⋯,I_{N,k}]^T$     (7)

and the component of matrices A,B and C are given, respectively, by

$A_{ll'}=1/(2l'+1) [l' δ_{l+1,l'}+(l'+1) δ_{l-1,l'}] $    (8a)

$B_{ll'}=(l' (l'+1))/(2l'+1) [δ_{l+1,l' }-(l'+1) δ_{l-1,l'}]$  (8b)

${{C}_{ll\text{ }\!\!'\!\!\text{ }}}={{\delta }_{ll\text{ }\!\!'\!\!\text{ }}}(1-\frac{(\omega {{a}_{l}})}{2l+1})$    (8c)

where δ is the Kronecker symbol; the components $D_l^k=δ_{0l}[1-ω(τ)] I_b (T)$ and $a_l$  is the scattering phase function coefficient. The Marshak boundary conditions are used in this study and written in matrix-vector form $A_0 F^1 (τ_0 )=D_0$ for inner and $A_L F^{N_L} (τ_L )=D_L$ for outer boundary [11, 18]

For transient heat transfer, the analytical solution of Eq. (7) for a single homogeneous and isothermal layer is similar to that developed by Ymeli and Kamdem [11] and Kamdem et al. [21] for planar media using double spherical harmonics method and discrete ordinates methods, respectively. The solution in each sub-layer is constructed by setting $F^k (τ)=R^k W^k (τ)$, where $R^k$ is the matrix of real eigenvector and $W^k (τ)$ is the vector of the characteristics intensity to be determined. Therefore, Eq. (7) is reduced to

$(dW^k)/dτ+Λ^k W^k=S^k=(AR^k )^{-1} D^k $    (9)

with $Λ^k$, the matrix of non-zero real eigenvalues  λ of the system. Each decoupled component of characteristics radiative intensity $W^k$ can be solved analytically and independently as [21]

$W^k (τ)=exp[λ^k (τ_B^k-τ)] C^k+\widetilde{\exp }[λ^k (τ_B^k-τ)] S^k$    (10)

with $τ_B^k=τ_l^k+((1-k_a )(τ_r^k-τ_l^k))⁄2$, where $τ_l^k$ and $τ_r^k$ are the optical depth of the left and right interfaces of the $k^{th}$ layer. The constant $k_a$ is $|λ^k |⁄λ^k$ and the diagonal matrices exp and $\widetilde{\exp }$ are defined as

$exp\left( {{\lambda }^{k}}(\tau _{B}^{k}-\tau ))=diag({{e}^{\lambda _{0}^{k}(\tau _{B}^{k}-\tau )}},...,{{e}^{\lambda _{N}^{k}(\tau _{B}^{k}-\tau )}} \right)$  (11)

$\widetilde{\exp }[λ^k (τ_B^k-τ)]=(Λ^k )^(-1).(I-exp{λ^k (τ_B^k-τ)})$  (12)

where I is the identity matrix with the same size with exp. The determination of the constant $C^k$ from boundary conditions and continuity condition is reduced to the linear system MX=N [11] where components are given by

${{M}_{(i,j)}}\text{=}\left\{ \begin{matrix}   \begin{matrix}   {{A}_{0}}{{R}^{1}}exp[{{x}_{1}}(\tau )] & i=j=1  \\\end{matrix}  \\   \begin{matrix}   {{R}^{j}}exp[{{x}_{j}}(\tau )]{{\delta }_{(i-1,j)}}-{{R}^{(j+1)}}exp[{{x}_{(j+1)}}(\tau ]{{\delta }_{(i,j)}} & {}  \\\end{matrix}  \\   \begin{matrix}   {{A}_{L}}{{R}^{({{N}_{L}})}}exp[{{x}_{({{N}_{L}})}}(\tau ) & i=j={{N}_{L}}  \\\end{matrix}  \\\end{matrix} \right.$    (13a)

${{N}_{i}}=\left\{ \begin{matrix}   {{D}_{0}}-{{A}_{0}}{{R}^{1}}\widetilde{exp}[{{x}_{1}}(\tau )]{{S}^{1}}j=1  \\   {{R}^{(j+1)}}\widetilde{exp}[{{x}_{(j+1)}}(\tau )]{{S}^{(j+1)}}-{{R}^{j}}\tau [{{x}_{j}}(\tau )]{{S}^{j}}  \\   {{D}_{L}}-{{A}_{L}}{{R}^{({{N}_{L}})}}\widetilde{exp}[{{x}_{({{N}_{L}})}}(\tau )]{{S}^{({{N}_{L}})}}j={{N}_{L}}  \\\end{matrix} \right.$ (13b)

where $X=[C^1,C^2,…,C^{N_L}]^T$ is the vector containing the integration constant of all the layers, $N_L$ is the total number of sub-layers and $x_i (τ)=λ^i (τ_B^i-τ_i )$. The linear system is solved with LU factorization and the final solution in each layer is given by

$F^k (τ)=R^k exp[λ^k (τ_B^k-τ)] C^k+R^k\widetilde{\exp }[λ^k (τ_B^k-τ)] S^k $    (14)

The developed methodology presents the advantages that the volumetric radiation can be computed directly at a particular selected location with reasonable computational time allowing the coupled conduction and radiation heat transfer. The radiative flux $q_R$ is obtained from $q_R^k (τ)=4πI_(1,k) (τ)/3$.

In this study, the temperature is initially guessed using the distribution at the previous time step and the assumption of multilayers structure of the medium where temperature and  are discretes and constants. This assumption has been also considered by Tan et al. [22] for radiation and Fourier’s conduction in graded index planar medium.

In addition, each layer has been splitting into nodes (including boundaries) to represent the solution given in Eq. (14). All these nodes are used to compute the solution of the conduction problem in the medium as a single layer.

3.2 Lattice Boltzmann method formulation (LBM)

The LBM is a computational method based on a mesoscopic description of heat flow, which better describes and captures sharp discontinuities of physical problems than the finite difference/volume methods [11, 23]. In order to simplify the energy equation with known radiative energy, we consider the dimensionless parameters as follows

$\eta =r/{{l}_{c}};\text{ }\xi ={{C}_{v}}t/{{l}_{c}},{{I}^{*}}=\frac{{{I}_{0}}(\eta )}{{{\sigma }_{B}}T_{h}^{4}/\pi },{{N}_{c}}=\frac{{{\kappa }_{e}}{{k}_{0}}}{4{{\kappa }_{B}}T_{h}^{3}}$   (15)

where $l_c=k_0⁄(ρC_p C_v )$ is the characteristic length of the problem, T_h is the reference temperature and N_c is the conductive-radiative parameter known also as Planck Number, the dimensionless time, distance, temperature $Θ=T/T_h$, and conductive heat flux $q_C^*=q_C/ρC_p T_h C_v$ are $ξ, η, Θ$, and $q_C^*$ respectively. Using this set of variables, Eqs. (1) and (2) are rewritten, as

$∂Θ/∂ξ+(∂q_c^*)/∂η=-a (q_C^*)/η-(κ_e l_c^2)/N_c  ∇^* q_R$  (16a)

$∂q_C^*/∂ξ+∂Θ/∂η=-q_C^*/u(Θ)$  (16b)

with $∇^* q_R  =κ_e [1-ω]{n^2 Θ^4-I^* }$ and $u(Θ)$ is a given function of temperature dependent thermal conductivity. The modified LBM equation due to radiation and non-Fourier is

$\frac{{{f}_{i}}}{\partial _{\xi }^{{}}}+{{\vec{e}}_{i}}.\nabla {{f}_{i}}=\frac{f{{_{i}^{e}}^{q}}-{{f}_{i}}}{{{\tau }_{M}}}-\frac{({{\kappa }_{e}}l_{c}^{2})}{2{{N}_{c}}}{{\nabla }^{*}}{{q}_{R}}-\frac{1}{2}\left( \frac{a}{\eta }/\hat{\eta }+\frac{{{e}_{i}}}{u(\theta )}) \right)q_{C}^{*}$   (17)

where $τ_M$ is the relaxation time, $f_i^{eq}$ is the equilibrium distribution function, $f_i$ is the particle distribution and $e ⃗_i$ is the moving velocity of the particles. For D1Q2 model, particles moves with two opposites unit velocities, $e_1  ⃗$  and $e_2 ⃗$. The relaxation time and vectors are $τ_M=u(Θ)+∆ξ/2$ [11, 23, 24]. After discretization, Eq. (17) is written as

${{f}_{i}}(\vec{\eta }+{{\vec{e}}_{i}}\Delta \xi ,\Delta +\Delta \xi )-{{f}_{i}}+\frac{\Delta \xi }{{{\tau }_{M}}}[f{{_{i}^{e}}^{q}}-{{f}_{i}}]+\frac{{{\kappa }_{e}}l_{c}^{2}}{2{{N}_{c}}}\Delta \xi {{\nabla }^{*}}{{q}_{R}}=\frac{\Delta \xi }{2}\left( \frac{a}{\eta }\hat{\eta }+\frac{{{{\vec{e}}}_{i}}}{u(\theta )} \right)q_{C}^{*}$    (18)

The temperature, conductive heat flux and equilibrium distributions are obtained from [11, 23, 24].

$Θ(η ⃗,ξ)=f_1 (η ⃗,ξ)+f_2 (η ⃗,ξ) $   (19a)

$q_C^* (η ⃗,ξ)=(e_1 ) ⃗.f_1 (η ⃗,ξ)+(e_2 ) ⃗.f_2 (η ⃗,ξ)$  (19b)

$f_i^{eq} (η ⃗,ξ)={Θ(η ⃗,ξ)+e ⃗_i.q_C^* (η ⃗,ξ)}/2$    (19c)

3.3 Quadrupole method

To solve the partial differential heat conduction equation, many different methods could be used but it is also very convenient to apply the Laplace transform, which transforms the partial differential equation into an ordinary differential equation in the Laplace domain. Let introduce p the Laplace variable and note To solve the partial differential heat conduction equation, many different methods could be used but it is also very convenient to apply the Laplace transform, which transforms the partial differential equation into an ordinary differential equation in the Laplace domain. Let introduce p the Laplace variable and note $L(f)=f ̅(p)=∫_0^{+∞}f(z,t)exp⁡(-pt)dt$.

$(d^2 Θ ̅)/(dz^2 )-τJ/k_C    (dΘ ̅)/dz-(ρC_p)/k_C  p(1+pτ_cv ) Θ ̅(z,p)+J^2/(pσk_C )=0$           (20)

The roots of the characteristic associate equation are

${{\gamma }_{1,2}}=\frac{\tau J}{2{{k}_{C}}}\pm \sqrt{{{(\frac{\tau J}{2{{k}_{C}}})}^{2}}+\frac{\rho {{C}_{p}}}{{{k}_{C}}}p(1+p{{\tau }_{cv}})}$   (21)

As a consequence, the solution of the Eq. (20) is

$Θ ̅(z,p)=ξ_1  exp⁡(γ_1 z)+ξ_2 exp(γ_2 z)+ϑ ̅(z)$  (22)

where $ϑ ̅(z)$ is a particular solution of Eq. (20) and $ξ_1, ξ_2$ are two constants depending on the boundary conditions. Two cases corresponding to different initial conditions are investigated to determine $ ϑ ̅(z)$:

IC0: the initial temperature within the leg is constant $T(z,t=0)=T_∞$ then it is obvious that

$\bar{\vartheta }(z)=\frac{{{J}^{2}}}{\rho {{C}_{p}}{{p}^{2}}\sigma (1+p{{\tau }_{c}}v)}+\frac{{{T}_{\infty }}}{p}$   (23)

IC1: the initial temperature is linear

$\bar{\vartheta }(z)=\frac{1}{p}(\frac{{{T}_{L}}-{{T}_{0}}}{L})z+\frac{\frac{{{J}^{2}}}{\sigma }+\frac{\tau J}{L}/({{T}_{L}}-{{T}_{0}})}{\rho {{C}_{p}}{{p}^{2}}(1+p{{\tau }_{c}}v)}+\frac{{{T}_{0}}}{p}$      (24)

Let for instance determine $ξ_1$ and $ξ_2$ in the most common case that is to say the thermoelectric thin layer is at the temperature $T_0$ and the side at $z=L$ is then at the temperature $T_L$. The equations corresponding to these boundary conditions are in the Laplace domain:

$Θ ̅(z=0,p)=0 ;   Θ ̅(z=L,p)=(T_L-T_0 )/p$  (25)

Then $ξ_1$ and $ξ_2$ are determined as

${{\xi }_{1}}=\frac{\frac{{{J}^{2}}}{\rho {{C}_{p}}{{p}^{2}}\sigma (1+p{{\tau }_{c}}_{v})}(1-exp({{\gamma }_{2}}L)-\frac{({{T}_{L}}-{{T}_{0}})}{p})}{exp({{\gamma }_{2}}L)-{{\xi }_{2}}exp({{\gamma }_{1}}L)}$    (26a)

${{\xi }_{2}}=\frac{\frac{{{J}^{2}}}{\rho {{C}_{p}}{{p}^{2}}\sigma (1+p{{\tau }_{c}}_{v})}(exp({{\gamma }_{2}}L)-1+\frac{({{T}_{L}}-{{T}_{0}})}{p})}{exp({{\gamma }_{2}}L)-{{\xi }_{2}}exp({{\gamma }_{1}}L)}$   (26b)

Now, let express the heat flux which is a linear combination of the temperature and the derivative of the temperature (where $α_s$ is the Seebeck coefficient)

$φ=α_s IΘ-k_C A ∂Θ/∂z$   (27)

Considering Eqs (22) and (27), the Laplace transform of the heat flux is:

$\bar{\varphi }\left( z \right)=\left( {{\alpha }_{s}}I-{{k}_{C}}A{{\gamma }_{1}} \right)\text{ }{{\xi }_{1}}~exp\left( {{\gamma }_{1}}z \right)+\left( {{\alpha }_{s}}I-{{k}_{C}}A{{\gamma }_{2}} \right)\text{ }{{\xi }_{2}}\text{ }exp\left( {{\gamma }_{2}}z \right)+\overline{{{\varphi }_{2}}}\left( z \right)$    (28)

$\overline{{{\varphi }_{2}}}\left( z \right)\text{=}~\left( {{\alpha }_{s}}I.\bar{\vartheta }(z)-{{k}_{C}}A\frac{d\bar{\vartheta }(z)}{dz} \right)$    (29)

From a mathematical point of view, the quadrupole method belongs to the class of analytical unified exact explicit method for solving linear partial differential equations in simple geometries. H.S. Carlaw [25] first presented this approach for the conduction of heat in solids. The Laplace heat flux and temperature are analogue of the electrical current and electrical potential. To summarize, this method provides a transfer matrix for the medium that linearly links the input temperature–heat flux column vector at one side and the output vector at the other side. Let now determine the quadrupole corresponding to the case studied. Thanks to Eqs. (22) and (28)

$\text{ }\left[ \begin{matrix}   \bar{\Theta }\left( {{z}_{i}} \right)  \\   \bar{\varphi }\left( {{z}_{i}} \right)  \\\end{matrix} \right]\text{=}{{M}_{p,i}}\left( \begin{matrix}   {{\xi }_{1}}  \\   \xi {}_{2}  \\\end{matrix} \right)+{{U}_{p,i}}$ (30)

where

111111.png

(31a)

$Mp.i=\left[ \begin{matrix}   \bar{\vartheta }({{z}_{i}})  \\   {{a}_{s}}I\bar{\vartheta }({{z}_{i}})-{{k}_{c}}A\bar{\vartheta }({{z}_{i}})/dz  \\\end{matrix} \right]$ (31b)

Considering the Eq. (26) at the both sides of the thermoelectric film, the quadrupole formulation of the problem is directly given by

$\text{ }\left[ \begin{matrix}   \bar{\Theta }\left( 0 \right)  \\   \bar{\varphi }\left( 0 \right)  \\\end{matrix} \right]\text{=}{{M}_{p,0}}{{M}_{p,L}}^{-1}\left[ \begin{matrix}   \bar{\Theta }\left( L \right)  \\   \bar{\varphi }\left( L \right)  \\\end{matrix} \right]-{{M}_{p,0}}{{M}_{p,L}}^{-1}{{U}_{p,L}}+{{U}_{p,0}}$ (32)

Thanks to this semi-analytical formulation, it is possible to plot the transient temperature along the thermoelectric film.

The transient coupled non-Fourier conduction with radiation heat transfer in an inhomogeneous, scattering and gray medium with graded index is investigated. The analytical layered solution of radiation is developed based on linear algebra approach in complex media while the lattice Boltzmann method is used to solve the non-linear hyperbolic energy equation and compared with benchmark solutions. The contribution of radiation on the transient temperature distributions in the medium was graphically examined using three radiation parameters: temperature dependent thermal conductivity, single scattering albedo and conduction-radiation parameter. The obtained results show good agreement with benchmark solutions of radiative heat flux, steady state and transient temperature distributions. It was found that the radiation effect is more pronounced for high value of refractive index or low values of the single scattering albedo; conduction-radiation parameter while this contribution is less effective for anisotropic scattering coefficient. Through these calculations, it was found that although the hyperbolic sharp wave front of non-Fourier conduction becomes smoother when strong radiative heat transfer is taken into account, it can be damped or emphasized with the effect of thermal conductivity and the scattering albedo. The present study demonstrates that the proposed analytical layered solution for radiation coupled to Lattice Boltzmann method for non-linear hyperbolic conduction is accurate and suitable for combined radiation-conduction problems in complex media. In the case of thermoelectric elements, the development of a specific quadrupole applied to thermoelectricity allows to obtain the transient behavior of a film under several boundary conditions, to take into account what happens at very short times and to investigate the effect of relaxation times on the thermogram.