## Abstract

Analytical modeling of thermal conduction in a multilayer body is of practical importance in several engineering applications such as microelectronics cooling, building insulation, and micro-electromechanical systems. A number of analytical methods have been used in past work to determine multilayer temperature distribution for various boundary conditions. However, there is a lack of work on solving the multilayer thermal conduction problem in the presence of spatially varying convective heat transfer boundary condition. This paper derives the steady-state temperature distribution in a multilayer body with spatially varying convective heat transfer coefficients on both ends of the body. Internal heat generation within each layer and thermal contact resistance between layers are both accounted for. The solution is presented in the form of an eigenfunction series, the coefficients of which are shown to be governed by a set of linear, algebraic equations that can be easily solved. Results are shown to be in good agreement with numerical simulation and with a standard solution for a special case. The model is used to analyze heat transfer for two specific problems of interest involving spatially varying convective heat transfer representative of jet impingement and laminar flow past a flat plate. In addition to enhancing the theoretical understanding of multilayer heat transfer, this work also contributes toward design and optimization of practical engineering systems comprising multilayer bodies.

## 1 Introduction

Thermal conduction in a multilayer body is relevant for the cooling and thermal management of several engineering systems [1,2]. For example, thermal management of packaged semiconductor chips typically involves heat transfer through a stack of multiple, thermally dissimilar materials. The boundary conditions at the two ends of a semiconductor chip package are also usually quite dis-similar. Typically, one end of the chip is attached to a board for electronic communication, resulting in poor rate of heat removal from that end. In contrast, the other end is typically used for heat spreading and removal using a heat sink. Another example of the practical importance of thermal conduction in multilayer bodies is in multilayer materials used for building insulation [3]. Finally, several microsystems comprise stacked thin films of multiple materials, where heat transfer through the stack plays a key role in determining the thermal characteristics of the microsystem [4].

Analytical modeling of heat transfer in such multilayer systems is clearly important for design and optimization, particularly when experiments-driven optimization may be expensive and time-consuming. The impact of key parameters such as thermal properties of various layers, heat generation rates and thermal contact resistances between layers on the temperature distribution is important to determine.

Significant literature exists on analytical modeling of thermal conduction in multilayer bodies. Solutions for multilayer problems in one-dimensional bodies have been presented [5,6]. Complex variables and residue theory methods have been used [7], but the associated computation is very tedious for more than a few layers [8]. Quasi-orthogonal methods have also been presented for a multilayer problem [9]. An exact solution for heat conduction in multilayer spherical composite laminates has been developed using separation of variables method [10]. The asymmetric heat conduction problem in a multilayer annulus with time-dependent boundary conditions has been solved using the finite integral transform method [11]. Using the Laplace transform method, semi-analytical solution for transient heat conduction in multilayer body within the framework of dual phase lag model has been presented [12]. Green’s function approach has been used to determine the temperature distribution in a 3D two-layer orthotropic slab [2]. The shifting function method and orthogonal expansion technique have been combined to determine temperature distributions of multilayered media with time-dependent heat transfer coefficient [13]. Analytical solution for transient heat conduction in a multilayer slab including various combinations of boundary conditions has been derived using a combination of superposition approach, separation of variables method, and orthogonal expansion technique [14]. Transient convective and radiative cooling over a multilayer composite slab has been modeled using the lumped parameter approach [15]. Finally, the thermal resistance network method has been used to model heat conduction problems in multilayer media [16,17]. Research has also been reported on modeling of the effects of thermal contact resistance in multilayer structures [1824]. In addition to prediction of the temperature distribution in such systems, several problems in this direction also relate to the use of inverse methods for estimation of time-dependent [20] and time- and space-dependent [23] thermal contact resistance. Both analytical [18,24] and numerical [19,21,22] approaches have been reported.

While boundary conditions of all three kinds have been analyzed in the literature summarized above, boundary condition of the third kind may be the most realistic, and therefore, of particular interest. Most of the literature in this direction is based on an assumption of a uniform constant convective heat transfer coefficient, h, applied over the entire boundary. Despite the simplification this assumption offers, a constant heat transfer coefficient may not be completely representative for several engineering systems. Spatial variation in h may occur due to spatially distributed fluid flow responsible for convective heat removal from the boundary. For example, in the case of thermal management by impinging jets, convective heat transfer on the surface has strong spatial distribution, being very large close to the site of impingement and reducing significantly farther outward [25]. Even laminar flow past a flat plate entails significant variation in h [26]. In light of this, accounting for spatial variation in convective heat transfer boundary conditions is important, but has not attracted sufficient research. The general solutions presented in the past [1,2] assume a constant convective heat transfer coefficient. A few mathematical treatments that account for spatial variation in h are available, but are limited only to single-layer or two-layer bodies [27,28]. Extension to a general multilayer body is desirable since several engineering and biomedical systems undergoing heat transfer due to such boundary conditions are multilayer in nature—a three-dimensional integrated circuit (3D IC) and skin are two common examples. For a homogeneous body, spatial variation in the convective boundary condition has been handled by considering a finite number of terms of an eigenfunction expansion, and then accounting for the spatial variation in the convective heat transfer coefficient by deriving a set of linear algebraic equations in the coefficients of the series [29,30]. This approach has been used for deriving analytical solutions for thermal conduction problems in a cylinder [29], fin [30], and sphere [31]. In addition, thermal conduction in a semiconductor chip with spatially varying convective boundary condition imposed on one end has also been addressed, for both homogeneous [32] and multilayer chip [33].

This paper presents a derivation of the temperature distribution in a multilayer body with distinct, spatially varying convective boundary conditions imposed on both ends of the body. This scenario models several engineering systems, and therefore, is of practical importance. The temperature distribution is written in the form of distinct series solutions for each layer based on eigenvalues from the side boundary conditions. Equations that relate coefficients for each layer with those of the bottom-most layer are derived. In turn, the coefficients for the bottom-most layer are determined by deriving a set of linear algebraic equations based on the two spatially varying convective heat transfer coefficients. The closed-form analytical solution represents an improvement over numerical simulations because the analytical solution does not require meshing or numerical simulation software. In addition, the analytical solution may be computationally more efficient, and may help provide physical insights into the nature of the heat transfer that is not possible through numerical simulations.

## 2 Analytical Modeling

Consider thermal conduction in an M-layer, composite body with spatially varying heat convective coefficients h1(x) and hM(x) on the two ends, as shown in Fig. 1. The thickness and thermal conductivity of each layer are (zi+1-zi) and ki, respectively. Uniform heat generation Qi occurs in each layer. Further, an interlayer thermal contact resistance Ri exists between the (i−1)th and ith layers. Each layer is assumed to be adiabatic at the two ends in the x direction. In a rectangular coordinate system, the governing steady-state energy conservation equation for the ith layer can be expressed as
$ki∂2Ti∂x2+ki∂2Ti∂z2+Qi=0(i=1,2,3,…,M)$
(1)
Fig. 1
Fig. 1
Close modal
Boundary conditions for this problem are
$∂Ti∂ x=0 at x=0,a$
(2a)
$k 1∂T1∂ z=h1(x)(T1−T∞,1) at z=0$
(2b)
$−kM∂TM∂ z=hM(x)(TM−T∞,M) at z=zM+1$
(2c)
Interfacial boundary conditions accounting for temperature and heat flux continuity are
$Ti−1(x,z)=Ti(x,z)−kiRi∂Ti∂z at z=zi$
(2d)
$ki−1∂Ti−1∂ z=ki∂Ti∂ z at z=zi (i=2,3,…,M)$
(2e)
In order to derive a solution for the steady-state temperature distribution in each layer, the following substitution is made in Eqs. (1)2(e):
$Ti(x,z)=θi(x,z)−Qi2kiz2$
(3)
This results in the following equations for θi(x,z):
$∂2θi∂x2+∂2θi∂z2=0$
(4)
Subject to
$∂θi∂ x=0 at x=0,a$
(5a)
$k 1∂θ1∂ z−h1(x)θ1(x,z)=F1(x) at z=0$
(5b)
$kM∂θM∂ z+hM(x)θM(x,z)=FM(x) at z=zM+1$
(5c)
$θi−1(x,z)=θi(x,z)−kiR i∂θi∂z+zi(Qi−12ki−1zi−Qi2kizi+R iQi) at z=zi$
(5d)
$ki−1∂θi−1∂ z+zi(Qi−Qi−1)=ki∂θi∂ z at z=zi$
(5e)
where
$F1(x)=−T∞,1h1(x)$
(6a)
$FM(x)=QMzM+1+QMzM+122kMhM(x)+T∞,MhM(x)$
(6b)
Since the boundary conditions in the x direction are adiabatic, the solution for$θi(x,z)$may be expressed in the form of the following Fourier cosine series:
$θi(x,z)=Ci ,0(z)+∑n=1∞Ci,n(z)cos(nπxa)$
(7)

Since Eq. (7) already satisfies Eq. (5a), the coefficients $Ci,0(z)$ and $Ci,n(z)$ must be determined in order to satisfy the governing Eq. (4) and boundary conditions (5b)(5e).

Substituting Eq. (7) in Eq. (4) results in
$C″i ,0(z)+∑n=1∞C″i ,n(z)cos(nπxa)−∑n=1∞Ci ,n(z)(nπa)2 cos(nπxa)=0$
(8)
Separately examining Ci,0 and Ci,n results in the following expressions:
$Ci, 0(z)=A i, 0+Bi ,0z$
(9a)
$Ci ,n(z)=A i ,n exp(nπza)+Bi ,n exp(−nπza)$
(9b)

where Ai,0, Bi,0, Ai,n, and Bi,n are unknown coefficients for each layer i =1,2,…,M and for each eigenvalue n =1,2,3…. These unknown coefficients are determined by requiring the expression for θi(x,z) given by Eq. (7) to satisfy the boundary and interfacial conditions (5b)(5e).

By substituting Eq. (7) into Eqs. (5d) and (5e), and using Eqs. (9a) and (9b), one may obtain
$A i,0=A 1,0+(1−k1k2)z2B 1,0+k1∑j=3i(1kj−1−1kj)zjB 1,0+k1∑j=2iRjB 1,0−∑j=2i(Qj2kj−Qj−1kj+Qj−12kj−1)zj2−∑j=2iRjQj−1zj−∑j=3i(1kj−1kj−1)zj∑k=2j−1(Qk−Qk−1)zk+∑j=3iRj∑k=2j−1(Qk−Qk−1)zk$
(10a)
$Bi,0=k1kiB1,0+1ki∑j=2i(Qj−Qj−1)zj$
(10b)
Further, $A i,n$ and $Bi,n$ can be expressed recursively as
$A i,n=12[1+(1+φi)ki−1ki]A i−1,n+12[1−(1+φi)ki−1ki]B i−1,n exp(−2nπzia)$
(11a)
$Bi,n=12[1+(1−φi)ki−1ki]B i−1,n+12[1−(1−φi)ki−1ki]A i−1,n exp(2nπzia)$
(11b)
where
$φi=nπkiR ia$
(12)

Equations (11a) and (11b) can be further simplified to write the coefficients for each layer, Ai,n and Bi,n in terms of coefficients for the bottom-most layer, $A 1,n$and $B 1,n$ as follows:

$A i,n=ηi,AA 1,n+μi,AB 1,n$
(13a)
$Bi,n=ηi,BA 1,n+μi,BB 1,n$
(13b)

In Eq. (13), ηi,A, ηi,B, μi,A and μi,B are given by the following coupled, recursive relationships:
$ηi,A=12[1+(1+φi)ki−1ki]ηi−1,A+12[1−(1+φi)ki−1ki]exp(−2nπzia)ηi−1,B$
(14a)
$μi,A=121+(1+φi)ki−1kiμi−1,A+12[1−(1+φi)ki−1ki]exp(−2nπzia)μi−1,B$
(14b)
$ηi,B=12[1−(1−φi)ki−1ki]exp(2nπzia)ηi−1,A+12[1+(1−φi)ki−1ki]ηi−1,B$
(14c)
$μi,B=12[1−(1−φi)ki−1ki]exp(2nπzia)μi−1,A+12[1+(1−φi)ki−1ki]μi−1,B$
(14d)

Initial values for these recursive relationships are given by $η1,A=1$, $μ 1,A=0$, $η1,B=0$, $μ 1,B=1$.

The coefficients for the Mth layer are of particular interest. Setting $i=M$ in Eqs. (10a), (10b), (13a), and (13b) yields
$A M,0=A 1,0+(1−k1k2)z2B 1,0+k1∑j=3M(1kj−1−1kj)zjB 1,0+k1∑j=2MRjB 1,0+f1$
(15a)
$BM,0=k1kMB1,0+f2$
(15b)
$AM,n=ηM,AA 1,n+μM,AB 1,n$
(16a)
$BM,n=ηM,BA 1,n+μM,BB 1,n$
(16b)
where
$f1=−∑j=2M(Qj2kj−Qj−1kj+Qj−12kj−1)zj2−∑j=2MRjQj−1zj−∑j=3M(1kj−1kj−1)zj∑k=2j−1(Qk−Qk−1)zk+∑j=3MRj∑k=2j−1(Qk−Qk−1)zk$
(17)
$f2=1kM∑j=2M(Qj−Qj−1)zj$
(18)

Equations (10a) and (10b) and (13a) and (13b) express all the coefficients needed to define the temperature solution – Ai,0, Bi,0, Ai,n and Bi,n—in terms of $A 1,0$, $B 1,0$ as well as $A 1,n$ and $B 1,n$.

In order to determine $A 1,0$, $B 1,0$ as well as $A 1,n$ and $B 1,n$, and thus complete the solution, Eq. (7) is separately substituted into the boundary conditions involving spatially varying heat transfer coefficients, Eqs. (5b) and (5c). Using equations (9a) and (9b), with i = M as well as (15b)(16b), one can obtain
$k1B 1,0+k1∑n=1∞nπa(A 1,n−B 1,n) cos(nπxa)−h1(x)A 1,0−h1(x)∑n=1∞(A 1,n+B 1,n) cos(nπxa)=F1(x)$
(19)
and
$k1B 1,0+kMf2+f1hM(x)+k1B 1,0hM(x)∑j=2MRj+k1kMB1,0zM+1hM(x)+f2zM+1hM(x)+kM∑n=1∞nπa[(ηM,AA 1,n+μM,AB 1,n)exp(nπzM+1a)−(ηM,BA 1,n+μM,BB 1,n)exp(−nπzM+1a)] cos(nπxa)+A 1,0hM(x)+(1−k1k2)z2B1,0hM(x)+k1B1,0hM(x)∑j=3M(1kj−1−1kj)zj+hM(x)∑n=1∞[(ηM,AA 1,n+μM,AB 1,n)exp(nπzM+1a)+(ηM,BA 1,n+μM,BB 1,n)exp(−nπzM+1a)] cos(nπxa)=FM(x)$
(20)
Equations (19) and (20) involve the unknown coefficients $A 1,0$, $B 1,0$ as well as $A 1,n$ and $B 1,n$ that must be determined. While the summations in Eqs. (19) and (20) involve, in principle, an infinite number of eigenvalues, these summations in Eqs. (19) and (20) are truncated at $n=N$ for computation. Integrating Eqs. (19) and (20) from x =0 to x = a separately and rearranging, one may obtain
$−A 1,0∫0ah1(x)dx+k1aB 1,0−∑n=1NA 1,n∫0ah1(x)cos(nπxa)dx−∑n=1NB 1,n ∫0ah1(x)cos(nπxa)dx=∫0aF1(x)dx$
(21)
$A 1,0∫0ahM(x)dx+k1aB 1,0+k1B 1,0∑j=2MRj∫0ahM(x)dx+k1kMB1,0zM+1∫0ahM(x)dx+(1−k1k2)z2B1,0∫0ahM(x)dx+k1B1,0∑j=3M(1kj−1−1kj)zj∫0ahM(x)dx+∑n=1N[ηM,A exp(nπzM+1a)+ηM,B exp(−nπzM+1a)] A 1,n∫0ahM(x)cos(nπxa)dx+∑n=1N[μM,A exp(nπzM+1a)+μM,B exp(−nπzM+1a)] B 1,n∫0ahM(x)cos(nπxa)dx=∫0aFM(x)dx−kMf2a−f1∫0ahM(x)dx−f2zM+1∫0ahM(x)dx$
(22)
Also, multiplying both sides of Eqs. (19) and (20) by $cos(iπxa)$ for $i=1,2,…,N$, and then integrating from x =0 to x = a and rearranging results in
$−A 1,0∫0ah1(x)cos(iπxa)dx+iπk12A 1,i−iπk12B 1,i−∑n=1NA 1,n∫0ah1(x)cos(nπxa)cos(iπxa)dx −∑n=1NB 1,n∫0ah1(x)cos(nπxa)cos(iπxa)dx=∫0aF1(x)cos(iπxa)dx$
(23)
$A 1,0∫0ahM(x)cos(iπxa)dx+k1B 1,0∑j=2MRj∫0ahM(x)cos(iπxa)dx+k1kMB1,0zM+1∫0ahM(x)cos(iπxa)dx+(1−k1k2)z2B1,0∫0ahM(x)cos(iπxa)dx+k1B1,0∑j=3M(1kj−1−1kj)zj∫0ahM(x)cos(iπxa)dx+iπkM2[ηM,A exp(iπzM+1a)−ηM,B exp(−iπzM+1a)]A 1,i+iπkM2[μM,A exp(iπzM+1a)−μM,B exp(−iπzM+1a)]B 1,i+∑n=1N[ηM,A exp(nπzM+1a)+ηM,B exp(−nπzM+1a)] A 1,n∫0ahM(x)cos(nπxa)cos(iπxa)dx+∑n=1N[μM,A exp(nπzM+1a)+μM,B exp(−nπzM+1a)] B 1,n∫0ahM(x)cos(nπxa)cos(iπxa)dx=∫0aFM(x)cos(iπxa)dx−f1∫0ahM(x)cos(iπxa)dx−f2zM+1∫0ahM(x)cos(iπxa)dx$
(24)

The (2 N +2) linear algebraic equations involving (2 N +2) unknowns, A1,0, B1,0, A1,i and B1,i (i =1,2,…,N) given by Eqs. (21)(24) can be easily solved using, for example, matrix inversion. This completes the solution, and the temperature distribution in the multilayer body is given by Eq. (3), along with Eqs. (7), (9a), (9b), (10a), (10b), (13a), and (13b).

## 3 Results and Discussion

Temperature distribution in a multilayer body is computed for a number of scenarios with spatially varying convective heat transfer coefficients imposed on both ends. While the model presented in Sec. 2 is capable of handling any general heat transfer coefficient, three specific heat transfer coefficient distributions of practical importance are considered. The first one is simply a step function where the value of h(x) is different in two or more different regions on the face. The second one represents convective heat transfer due to laminar flow past a flat plate, in which case, the local convective heat transfer coefficient h(x) is known to scale as x−1/2 [26]. The third one corresponds to jet impingement cooling, for which h(x) is known to vary strongly with x, and has been modeled by the following equation [25]:
$h(x)=hmax[1−Ptanh(γ(|x−x0|d−1.5))1+P]$
(25)

where x0 and d are the location and width of the jet, respectively. $P=hmax−hminhmax+hmin$, where hmax and hmin represent the maximum and minimum convective heat transfer coefficients due to the jet near the impingement site and away from it, respectively. γ models the sharpness of transition between hmin and hmax.

### 3.1 Effect of Number of Eigenvalues.

Since the solution for temperature distribution is derived in the form of an infinite series, a minimum number of terms must be determined in order to compute the solution. This involves a key tradeoff because increasing the number of terms may improve accuracy but also increase computational cost. This tradeoff is analyzed in detail for a representative, two-layer problem with spatially varying convective heat transfer on both ends. Jet impingement cooling represented by Eq. (25) is implemented on the bottom face. Values of parameters appearing in Eq. (25) are x0 = 15 mm, d =5 mm, hmin = 30 W/(m2K), hmax = 200 W/(m2K) and γ = 2.0. On the top face, convective heat transfer due to laminar flow past a flat plate is modeled, such that h(x) = C⋅x−1/2, where C =6.0 Wm−1.5K−1. The values of T∞,1 and T∞,2 are 50 K and 100 K, respectively. Each layer is 5 mm thick and 30 mm wide. Thermal conductivities of the two layers are 1 W/(mK) and 2 W/(mK), respectively. No internal heat generation is modeled, but an interlayer thermal contact resistance of 10−5 Km2/W is modeled.

Under these conditions, Figs. 2(a) and 2(b) plot computed temperature distributions on the bottom and top faces, respectively. In each case, the computed distribution for five different total number of terms is considered. Both plots show that the computed distribution is quite inaccurate when only one eigenvalue is considered. However, there is rapid convergence with increasing number of eigenvalues, and beyond five eigenvalues, there is very little change in the computed temperature distribution by considering additional terms. Temperature computed with 20 and 30 terms is nearly identical to each other, with less than 0.5% deviation. This shows that around 20 eigenvalues is sufficient for accurate temperature computation using the results presented in Sec. 2. Data in all subsequent Figures in this section are computed with 20 eigenvalues.

Fig. 2
Fig. 2
Close modal

### 3.2 Validation of Analytical Model.

In order to validate the theoretical model presented in Sec. 2, results are compared with numerical simulations for a representative three-layer problem. The numerical simulation is carried out in ansys-cfx. The three-layer geometry is modeled and meshed with 60,000 elements. Mesh sensitivity study is carried out to ensure that the results are independent of any further refinement. Solver control for this steady-state thermal conduction problem is set to either 500 iterations or a residual of 10−8. The simulation is found to reach the target residual after around 250 iterations. For this comparison, step change h(x) is assumed on both faces. On the bottom face, h1(x)=100 W/(m2K) when 0≤x < a/3, and h1(x)=200 W/(m2K) when x  a/3. On the top face, hM(x) = 50 W/(m2K) when 0≤x <2a/3, and hM(x) = 150 W/(m2K) when x 2a/3. The values of T∞,1 and T∞,3 are 50 K and 100 K, respectively. Each layer is 5 mm thick, with a width of a =20 mm. Thermal conductivity values for the three layers are 2 W/(mK), 4 W/(mK) and 1 W/(mK), respectively. Internal heat generation of 10,000 W/m3 is modeled in each layer, and an interlayer thermal contact resistance of 10−5 Km2/W is assumed. For these conditions, Fig. 3(a) plots the entire temperature distribution in the multilayer body and compares results from the analytical model presented in this work with numerical simulations. Figure 3(b) presents this comparison specifically for the top and bottom face temperatures as functions of x. These figures show very good agreement between the analytical model and numerical simulations. The maximum deviation between the two is less than 0.1%, thereby indicating the accuracy of the analytical model.

Fig. 3
Fig. 3
Close modal
For further validation of the theoretical model, a special case of constant convective heat transfer coefficient is considered, and computed temperature distribution is compared with a well-known solution for the special case. For a two-layer body with no internal heat generation, a given interlayer thermal contact resistance R, and constant heat transfer coefficients h1 and hM at the two ends, the temperature distribution is clearly one-dimensional, and is given by
$T1(z)=T∞,1+ΔTfz+k1ΔTh1f$
(26)
$T2(z)=T∞,M+k1ΔTk2fz−k1(1hM+z3k2)ΔTf$
(27)
where
$ΔT=T∞,M−T∞,1$
(28)
$f=z2+k1h1−k1k2z2+k1hM+k1k2z3+k1R2$
(29)

Figure 4 compares the temperature distribution determined using the analytical model with the well-known solution given by Eqs. (26) and (27). In this case, the two layers are both 5 mm thick, with thermal conductivities of 1 W/(mK) and 2 W/(mK), respectively. The value of R is 2 × 10−3 Km2/W, and the heat transfer coefficients are 50 W/(m2K) and 100 W/(m2K), respectively. The values of T∞,1 and T∞,2 are 100 K and 50 K, respectively. Based on these parameters, Fig. 4 shows excellent agreement between the analytical model and the well-known solution for this special case. The analytical model accurately captures the temperature gradient within each layer, as well as the temperature discontinuity at the interface due to the nonzero thermal contact resistance.

Fig. 4
Fig. 4
Close modal

### 3.3 Applications of the Analytical Model.

Applications of the analytical model for solving two practical heat transfer problems in multilayer bodies are presented next.

In the first problem, cooling of a three-layer heat generating body by jet impingement on both sides is considered. The body is 30 mm wide, each layer is 5 mm thick. Thermal conductivity values of 2 W/(mK), 5 W/(mK), and 1 W/(mK), respectively, for the three layers. Internal heat generation of 800,000 W/m3 is assumed, along with interfacial thermal contact resistance of 10−5 Km2/W. Jet parameters for the bottom jet are x0=15 mm, d =5 mm, hmin = 100 W/(m2K), hmax = 1000 W/(m2K), and γ = 2.0, whereas the values for the top jet are x0 = 25 mm, d =2 mm, hmin = 100 W/(m2K), hmax = 500 W/(m2K), and γ = 2.0. As evident from these values, the top jet is somewhat weaker and narrower, and located near the right end of the body, while the bottom jet is stronger, broader, and located at the center of the body. The ambient temperature for convective heat transfer is assumed to be 0 K for both jets so that the model effectively computes temperature rise in the body above ambient. For these parameters, Fig. 5(a) plots the temperature colormap in the body, and Fig. 5(b) plots temperature distribution along the top and bottom plates. These results are consistent with the expected temperature distributions based on the values of the parameters. For example, temperature along the bottom face is nearly symmetric due to the center location of the jet on the bottom face. On the other hand, the top face is coolest around x =25 mm, which is where the jet on the top face is centered. The widths of the temperature distributions for the top and bottom plates shown in Fig. 5(b) are consistent with the assumed values of d for the two jets. Within the body, the peak temperature occurs near the left end, which is farthest away from the two jets, and is close to the interface between the second and third layers.

Fig. 5
Fig. 5
Close modal

The second problem considers the heating of a three-layer body with laminar flow of a hot fluid on both faces of the body. The body is 30 mm wide, and each layer is 5 mm thick, with a thermal contact resistance of 10−5 Km2/W at each interface, no internal heat generation and thermal conductivities of 1 W/(mK), 2 W/(mK) and 5 W/(mK), respectively. The convective heat transfer coefficients on the two faces are taken to be h1(x) = C1x−1/2 and hM(x) = CMx−1/2, where C1 = 7.4  Wm−1.5K−1, and CM = 6.0  Wm−1.5K−1. These values roughly correspond to air flow with freestream velocities of 15 m/s and 10 m/s, respectively. The freestream temperatures for the two flows are taken to be 100 K and 50 K, respectively. Under such conditions, it is expected that there will be greater heat transfer along the bottom plate than the top face, and that, within each face, heat transfer will be greatest at small values of x. As the thermal boundary layer grows and the value of h decreases, the rate of heat exchange between the body and the fluid is expected to go down. The temperature distribution in the three-layer body is computed using the analytical model for the parameters listed above. Results are shown in Fig. 6, where Fig. 6(a) shows a colormap of the entire temperature distribution, while Fig. 6(b) plots temperature distribution as a function of x on the top and bottom faces. Results are along expected lines. The bottom face is hotter than the top face, as expected, due to the greater freestream temperature associated with convection. Due to the x−1/2 dependence of the convective heat transfer coefficient, Fig. 6(b) shows that temperature on the bottom face reduces with increasing x, which causes the bottom face temperature to go farther away from the T∞,1 = 100 K value. Similarly, the rate of convective heat transfer from the top face to the second fluid at T∞,M = 50 K reduces as x increases, and therefore, the top face temperature increases with x away from the freestream temperature of the second fluid.

Fig. 6
Fig. 6
Close modal

## 4 Conclusions

The capability of the analytical model presented here to account for spatially dependent convective heat transfer on both ends of the body facilitates rapid computation of temperature distribution in the multilayer body. The solution is presented in a simple series form, the coefficients of which are easy to compute. The analysis related to the number of eigenvalues needed helps balance accuracy with computational cost. In addition to improving our fundamental understanding of heat transfer in a multilayer body, the theoretical model presented in this paper may be helpful in analyzing and optimizing multilayer heat transfer in a number of engineering systems. In addition, the model may also help understand bioheat transfer in multilayer systems such as skin [34].

## Funding Data

• Key Project of Science of the Education Bureau of Henan Province (Grant No. 19B460005).

• Special Project of Basic Scientific Research Operating Expenses of Henan Polytechnic University (Grant No. NSFRF180427).

• China Scholarship Council (Funder ID: 10.13039/501100004543).

• National Science Foundation (Grant No. 1554183; Funder ID: 10.13039/100000001).

## Nomenclature

h =

convective heat transfer coefficient (W/(m2K))

k =

thermal conductivity (W/(mK))

M =

number of layers

N =

number of eigenvalues

Q =

internal heat generation rate (W/m3)

R =

thermal contact resistance (Km2/W)

T =

temperature (K)

x,z =

spatial coordinates (m)

## References

1.
Choobineh
,
L.
, and
Jain
,
A.
,
2015
, “
An Explicit Analytical Model for Rapid Computation of Temperature Field in a Three-Dimensional Integrated Circuit
, (3D IC),”
Int. J. Therm. Sci.
,
87
, pp.
103
109
.10.1016/j.ijthermalsci.2014.08.012
2.
Haji-Sheikh
,
A.
, and
Beck
,
J.
,
2002
, “
Temperature Solution in Multi-Dimensional Multi-Layer Bodies
,”
Int. J. Heat Mass Transfer
,
45
(
9
), pp.
1865
1877
.10.1016/S0017-9310(01)00279-4
3.
Kang
,
D.-H.
,
Lorente
,
S.
, and
Bejan
,
A.
,
2013
, “
Constructal Distribution of Multi‐Layer Insulation
,”
Int. J. Energy Res.
,
37
(
2
), pp.
153
160
.10.1002/er.1895
4.
Hui
,
P.
, and
Tan
,
H. S.
,
1994
, “
A Transmission-Line Theory for Heat Conduction in Multilayer Thin Films
,”
IEEE Trans. Compon., Packag. Manuf. Technol. Part B
,
17
(
3
), pp.
426
434
.10.1109/96.311793
5.
Hahn
,
D. W.
, and
Özişik
,
M. N.
,
2012
,
Heat Conduction
,
Wiley, New York
.
6.
Mikhailov
,
M. D.
, and
Özişik
,
M. N.
,
1994
,
Unified Analysis and Solutions of Heat and Mass Diffusion
,
Dover Publications
,
New York
.
7.
Carslaw
,
H. S.
, and
Jaeger
,
J. C.
,
1986
,
Conduction of Heat in Solids
,
Oxford Science Publications
,
London
.
8.
Mulholland
,
G. P.
, and
Cobble
,
M. H.
,
1972
, “
Diffusion Through Composite Media
,”
Int. J. Heat Mass Transfer
,
15
(
1
), pp.
147
160
.10.1016/0017-9310(72)90172-X
9.
Tittle
,
C. W.
,
1965
, “
Boundary Value Problems in Composite Media: Quasi‐Orthogonal Functions
,”
J. Appl. Phys.
,
36
(
4
), pp.
1486
1488
.10.1063/1.1714335
10.
Norouzi
,
M.
,
Amiri Delouei
,
A.
, and
Seilsepour
,
M.
,
2013
, “
A General Exact Solution for Heat Conduction in Multilayer Spherical Composite Laminates
,”
Compos. Struct.
,
106
, pp.
288
295
.10.1016/j.compstruct.2013.06.005
11.
Singh
,
S.
, and
Jain
,
P. K.
,
Rizwan-Uddin
,
2011
, “
Finite Integral Transform Method to Solve Asymmetric Heat Conduction in a Multilayer Annulus With Time-Dependent Boundary Conditions
,”
Nucl. Eng. Des.
,
241
, pp.
144
154
. 10.1016/j.nucengdes.2010.10.010
12.
,
K.
,
2009
, “
Semi-Analytical Solutions for the Dual Phase Lag Heat Conduction in Multilayered Media
,”
Int. J. Therm. Sci.
,
48
(
1
), pp.
14
25
.10.1016/j.ijthermalsci.2008.03.004
13.
Chiba
,
R.
,
2018
, “
An Analytical Solution for Transient Heat Conduction in a Composite Slab With Time-Dependent Heat Transfer Coefficient
,”
Math. Probl. Eng.
,
2018
, pp.
1
11
.10.1155/2018/4707860
14.
Wu
,
X.
,
Shi
,
J-y.
,
Lei
,
K.
,
Li
,
Y-P.
, and
Leslie
,
O.
,
2019
, “
Analytical Solutions of Transient Heat Conduction in Multilayered Slabs and Application to Thermal Analysis of Landfills
,”
J. Central South Univ.
,
26
(
11
), pp.
3175
3187
.10.1007/s11771-019-4244-y
15.
An
,
C.
, and
Su
,
J.
,
2011
, “
Improved Lumped Models for Transient Combined Convective and Radiative Cooling of Multi-Layer Composite Slabs
,”
Appl. Therm. Eng.
,
31
(
14–15
), pp.
2508
2517
.10.1016/j.applthermaleng.2011.04.016
16.
Noh
,
J.-H.
,
Kim
,
W.-G.
,
Cha
,
K.-U.
, and
Yook
,
S.-J.
,
2015
, “
Inverse Heat Transfer Analysis of Multi-Layered Tube Using Thermal Resistance Network and Kalman Filter
,”
Int. J. Heat Mass Transfer
,
89
, pp.
1016
1023
.10.1016/j.ijheatmasstransfer.2015.06.009
17.
Noh
,
J.-H.
,
Cha
,
K.-U.
,
Ahn
,
S.-T.
, and
Yook
,
S.-J.
,
2016
, “
Prediction of Time-Varying Heat Flux Along a Hollow Cylindrical Tube Wall Using Recursive Input Estimation Algorithm and Thermal Resistance Network Method
,”
Int. J. Heat Mass Transfer
,
97
, pp.
232
241
.10.1016/j.ijheatmasstransfer.2016.02.011
18.
Belghazi
,
H.
,
El Ganaoui
,
M.
, and
Labbe
,
J. C.
,
2010
, “
Analytical Solution of Unsteady Heat Conduction in a Two-Layered Material in Imperfect Contact Subjected to a Moving Heat Source
,”
Int. J. Therm. Sci.
,
49
(
2
), pp.
311
318
.10.1016/j.ijthermalsci.2009.06.006
19.
Frekers
,
Y.
,
Helmig
,
T.
,
Burghold
,
E. M.
, and
Kneer
,
R.
,
2017
, “
A Numerical Approach for Investigating Thermal Contact Conductance
,”
Int. J. Therm. Sci.
,
121
, pp.
45
54
.10.1016/j.ijthermalsci.2017.06.026
20.
Yang
,
Y.-C.
,
Lee
,
H.-L.
,
Chen
,
W.-L.
, and
Salazar
,
J. L.
,
2011
, “
Estimation of Thermal Contact Resistance and Temperature Distributions in the Pad/Disc Tribosystem
,”
Int. Commun. Heat Mass Transfer
,
38
(
3
), pp.
298
303
.10.1016/j.icheatmasstransfer.2010.11.005
21.
Sakkaki
,
M.
,
,
F.
,
Vajdi
,
M.
,
Pishgar
,
F.
,
Shokouhimehr
,
M.
, and
Shahedi Asl
,
M.
,
2019
, “
The Effect of Thermal Contact Resistance on the Temperature Distribution in a WC Made Cutting Tool
,”
Ceram. Int.
,
45
(
17
), pp.
22196
22202
.10.1016/j.ceramint.2019.07.241
22.
Shittu
,
S.
,
Li
,
G.
,
Zhao
,
X.
,
Ma
,
X.
,
Akhlaghi
,
Y. G.
, and
Fan
,
Y.
,
2020
, “
Comprehensive Study and Optimization of Concentrated Photovoltaic-Thermoelectric Considering All Contact Resistances
,”
Energy Convers. Manage.
,
205
, p.
112422
10.1016/j.enconman.2019.112422.
23.
Gaspar
,
J.
,
Rigollet
,
F.
,
Gardarein
,
J.-L.
,
Le Niliot
,
C.
, and
Corre
,
Y.
,
2014
, “
Identification of Space and Time Varying Thermal Resistance: Numerical Feasibility for Plasma Facing Materials
,”
Inverse Prob. Sci. Eng.
,
22
(
2
), pp.
213
231
.10.1080/17415977.2012.757314
24.
Fakoor-Pakdaman
,
M.
,
,
M.
,
Bagheri
,
F.
, and
Bahrami
,
M.
,
2015
, “
Optimal Time-Varying Heat Transfer in Multilayered Packages With Arbitrary Heat Generations and Contact Resistance
,”
ASME J. Heat Transfer
,
137
(
8
), p.
081401
.10.1115/1.4028243
25.
Sarkar
,
D.
,
Jain
,
A.
,
Goldstein
,
R.
, and
Srinivasan
,
V.
,
2016
, “
Corrections for Lateral Conduction Error in Steady-State Heat Transfer Measurements
,”
Int. J. Therm. Sci.
,
109
, pp.
413
423
.10.1016/j.ijthermalsci.2016.05.031
26.
Bergman
,
T. L.
,
Dewitt
,
D. P.
,
Lavine
,
A. S.
, and
Incropera
,
F. P.
,
2011
,
Fundamentals of Heat and Mass Transfer
,
Wiley
,
New York
.
27.
Özişik
,
M. N.
, and
Murray
,
R. L.
,
1974
, “
On the Solution of Linear Diffusion Problems With Variable Boundary Condition Parameters
,”
ASME J. Heat Transfer
,
96
(
1
), pp.
48
51
.10.1115/1.3450139
28.
Holy
,
Z. J.
,
1972
, “
Temperature and Stresses in Reactor Fuel Elements Due to Time- and Space-Dependent Heat-Transfer Coefficients
,”
Nucl. Eng. Des.
,
18
(
1
), pp.
145
197
.10.1016/0029-5493(72)90041-6
29.
Sarkar
,
D.
,
Shah
,
K.
,
Haji-Sheikh
,
A.
, and
Jain
,
A.
,
2014
, “
Analytical Modeling of Temperature Distribution in an Anisotropic Cylinder With Circumferentially-Varying Convective Heat Transfer
,”
Int. J. Heat Mass Transfer
,
79
, pp.
1027
1033
.10.1016/j.ijheatmasstransfer.2014.08.060
30.
Ma
,
S. W.
,
Behbahani
,
A. I.
, and
Tsuei
,
Y. G.
,
1991
, “
Two-Dimensional Rectangular Fin With Variable Heat Transfer Coefficient
,”
Int. J. Heat Mass Transfer
,
34
(
1
), pp.
79
85
.10.1016/0017-9310(91)90175-E
31.
Sarkar
,
D.
,
Haji-Sheikh
,
A.
, and
Jain
,
A.
,
2016
, “
Thermal Conduction in an Orthotropic Sphere With Circumferentially Varying Convection Heat Transfer
,”
Int. J. Heat Mass Transfer
,
96
, pp.
406
412
.10.1016/j.ijheatmasstransfer.2016.01.027
32.
Luhar
,
S.
,
Sarkar
,
D.
, and
Jain
,
A.
,
2017
, “
Steady State and Transient Analytical Modeling of Non-Uniform Convective Cooling of a Microprocessor Chip Due to Jet Impingement
,”
Int. J. Heat Mass Transfer
,
110
, pp.
768
777
.10.1016/j.ijheatmasstransfer.2017.03.064
33.
Zhou
,
L.
,
Parhizi
,
M.
, and
Jain
,
A.
,
2020
, “
Theoretical Modeling of Heat Transfer in a Multilayer Body With Spatially-Varying Convective Heat Transfer Boundary Condition
,”
Int. J. Heat Mass Transfer
, in review.
34.
Bhargava
,
A.
,
Chanmugam
,
A.
, and
Herman
,
C.
,
2014
, “
Heat Transfer Model for Deep Tissue Injury: A Step Towards an Early Thermographic Diagnostic Capability
,”
Diagn. Pathol.
,
9
(
1
), pp.
36
18
.10.1186/1746-1596-9-36