Analytical Solution for Temperature Distribution in a Multilayer Body With Spatially Varying Convective Heat Transfer Boundary Conditions on Both Ends

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 [18–24]. 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.

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).

where A_i,0, B_i,0, A_i,n, and B_i,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).

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

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

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

The (2 N + 2) linear algebraic equations involving (2 N + 2) unknowns, A_1,0, B_1,0, A_1,i and B_1,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).

where x₀ and d are the location and width of the jet, respectively. $P=hmax−hminhmax+hmin$⁠, where h_max and h_min 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 h_min and h_max.

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 x₀ = 15 mm, d = 5 mm, h_min = 30 W/(m²K), h_max = 200 W/(m²K) 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⁻¹. 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⁻⁵Km²/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.

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⁻⁸. 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, h₁(x)=100 W/(m²K) when 0≤x < a/3, and h₁(x)=200 W/(m²K) when x ≥ a/3. On the top face, h_M(x) = 50 W/(m²K) when 0≤x < 2a/3, and h_M(x) = 150 W/(m²K) 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/m³ is modeled in each layer, and an interlayer thermal contact resistance of 10⁻⁵Km²/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.

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⁻³Km²/W, and the heat transfer coefficients are 50 W/(m²K) and 100 W/(m²K), 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.

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/m³ is assumed, along with interfacial thermal contact resistance of 10⁻⁵Km²/W. Jet parameters for the bottom jet are x₀=15 mm, d = 5 mm, h_min = 100 W/(m²K), h_max = 1000 W/(m²K), and γ = 2.0, whereas the values for the top jet are x₀ = 25 mm, d = 2 mm, h_min = 100 W/(m²K), h_max = 500 W/(m²K), 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.

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⁻⁵Km²/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 h₁(x) = C₁⋅x^−1/2 and h_M(x) = C_M⋅x^−1/2, where C₁ = 7.4  Wm^−1.5K^−1, and C_M = 6.0  Wm^−1.5K⁻¹. 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.

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].

• 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).

h =

convective heat transfer coefficient (W/(m²K))

k =

thermal conductivity (W/(mK))

M =

number of layers

N =

number of eigenvalues

Q =

internal heat generation rate (W/m³)

R =

thermal contact resistance (Km²/W)

T =

temperature (K)

x,z =

spatial coordinates (m)