1. Introduction
Two or three-layer structures are used in a number of applications in engineering, such as in the field of thermoplastic shaping of polymers [
1] and in the estimation of the thermal properties of solid materials [
2,
3]. In the latter field, for example, an experimental apparatus consisting of a thin layer heater located between two specimens of the same material and thickness is employed in order to obtain measured temperature values. Additionally, this experimental configuration is usually reduced to a simplified two-layer configuration in which a surface layer, representing the thin heater, is in contact with a finite plate (specimen).
The complex thermal model of such structures is described by a set of energy equations coupled at the interface between each layer. Additionally, this thermal model may be simplified if the thermal diffusion within the surface layer (solid or fluid) is neglected, by introducing the assumption of lumped body. In particular, the boundary conditions appearing in the simplified model are modified so as to include different thermal effects due to the thin layer. Several “modified” boundary conditions exist in the heat conduction literature. One of these is the fourth kind (Carslaw) boundary condition [
4] which takes the thermal capacity of the thin layer into account; in fact, it derives from the application of the first law of thermodynamics to the thin layer. More general boundary conditions are represented by the fifth kind (Jaeger) boundary condition [
4] which not only takes the thermal capacity of the thin layer into account, but it also permits a convective heat flux from the thin layer to the surrounding ambient; and by the sixth kind boundary condition [
5] which is imposed when an imperfect thermal contact between the thin layer and an adjacent domain occurs.
Although many of multi-layer transient heat conduction problems have been studied in literature [
6,
7,
8,
9,
10], in several cases the derivations of these solutions are not provided as occurs, for example, for many problems treated in Ref. [
6]; also, for fewer cases the solution is not available at all.
In the present work one-dimensional multi-layer heat conduction problems involving boundary condition of high kinds (fourth and sixth kind) are analyzed. In order to easily identify which conductive problem is analyzed, the numbering system devised in [
11] is used. According to this notation the addressed problems may be listed as: (1) X40B1T00 for the case of a thin layer in perfect thermal contact with a semi-infinite body; (2) X60B1T00 when the contact between them is imperfect; (3) X42B10T00 for the case in which the thin layer is in perfect contact with a slab, insulated at the other boundary; and (4) X62B10T00 when an imperfect contact between them is taken into account. In particular, the notation X40B1T00 denotes a one dimensional transient heat conduction problem concerning a rectangular (by the “X”) semi-finite body (by the “0” in the “X40”) in perfect contact with a thin layer at the surface
x = 0 (fourth kind boundary condition by the “4” in “X40”) where a jump in heat flux is applied (by the B1); also, T00 stands for a zero initial temperature for both layers. Moreover, for the case of imperfect contact the “4” in “X40” is replaced by “6” (sixth kind boundary condition at the surface
x = 0). Furthermore, for the problem X42B10T00 (or X62B10T00 for imperfect contact) concerning a thin layer adjacent to a slab which is insulated at the other boundary, the “2” in “X42” (or in “X62”) denotes a heat flux boundary condition at
x =
L which is homogeneous by the “0” in “B10”.
Moreover, as the temperature solution to the X60B1T00 case is not available in the heat conduction literature, it is derived by means of the Laplace Transform (LT) method. Unlike this last case, the exact analytical solutions to the X42B10T00 and X62B10T00 problems are available in the literature, but their derivations are not given. Therefore, they are proposed by using the Green’s Function Solution Equation (GFSE) method, for the X42B10T00 case, and the LT method for the X62B10T00 case. Additionally, for both problems the temperature computational solutions are defined for short and large-times, as well as for the a quasi-steady state.
In addition, three-layer structures involving a thin layer with internal heat generation, located between two semi-infinite bodies or two finite plates (slabs for short) are discussed.
2. Thin Layer in Perfect Contact with a Semi-Infinite Body
Consider a semi-infinite rectangular one-dimensional (1D) body, initially at uniform temperature
Tin, in perfect thermal contact with a high-conductivity surface layer at the boundary
x = 0, as depicted in
Figure 1. At
a step change in heat flux is applied to the surface layer at its boundary
, which is at the same initial temperature
Tin. In particular, the layer is considered to be a lumped body; or, in other words, it is thin enough to be able to neglect the thermal gradients developing within it. Additionally, the thermal properties
k,
C and
Cf are considered temperature-independent. The problem notation is X40B1T11.
The governing equations for the addressed problem, in dimensionless form, are expressed as
where Equation (1b,c) represent the boundary condition of the fourth kind. In particular, Equation (1b) is obtained applying the first law of thermodynamics to the thin layer; while Equation (1c) indicates a perfect contact between the semi-infinite body and the thin layer (temperature continuity).
The dimensionless variables appearing in Equation (1) are defined as
where
L is a reference length, while
is the thickness of the thin layer. Additionally,
P denotes the heat capacity ratio. In dimensionless form, the problem notation becomes X40B1T00.
A well established exact analytical solution is available in the heat conduction literature ([
6], p. 306, Equation (12)) for the current problem. Additionally, this solution is discussed and derived in [
12], where a verified computer code is made available. It is
where
is the so-called Amos function, defined as
The is the scaled complementary error function. It is recommended for large values of , and as difficulties might arise during the evaluation of Equation (3b) in its former expression because of the positive argument of the exponential term. On the contrary, the scaled complementary error function (appearing in the latter expression) approaches zero for large values of , and .
2.1. Thin Layer between Two Semi-Infinite Bodies
The solution given by Equation (3) can also be used for the three-layer configuration shown in
Figure 2. In this case, a thin layer (still considered as a lumped body with temperature-independent properties) is in perfect contact with two semi-infinite bodies. The thin layer and the semi-infinite bodies are initially at uniform temperature
Tin, but at
t = 0 only the former is subject to a uniform and constant internal heat generation per unit of volume
.
In such a case, the governing equations in dimensionless form are
where Equation (4g) is obtained applying an energy balance to the thin layer; while Equation (4h,i) denote a perfect contact between the thin layer and the two semi-infinite bodies.
The dimensionless variables which appear in Equation (4) are defined as
For the sake of thermal symmetry, the three-layer structure depicted in
Figure 2 reduces to the two-layer configuration shown in
Figure 1 for perfect contact. In fact, by considering only an half-thickness of the thin layer (say
), and by noting that
in Equation (5) plays the same role of the surface heat flux
in Equation (2), the problem defined by Equation (4) reduces to the X40B1T00 case defined by Equation (1). Hence, the thermal field
within the semi-infinite body on the right-side of the three-layer configuration in
Figure 2 is still given by Equation (3a), which can also be used to describe the thermal field
on the left-side by simply replacing the variable
with
. As far as the thin film is concerned, its temperature is defined by Equation (3a) for
.
4. Thin Layer in Perfect Contact with a Slab
Consider now a slab with temperature-independent properties, initially at uniform temperature
Tin, in perfect contact with a thin layer at the boundary
x = 0, and insulated at the back side as depicted in
Figure 3. At
a step change in heat flux is applied to the thin film at
, which is at the same initial temperature
Tin. The problem notation is X42B10T11.
The mathematical formulation of the conductive problem, in dimensionless form, is reported below.
The dimensionless variable appearing in Equation (18) are defined in Equation (2), where
L is now the slab thickness. In dimensionless form, the problem notation becomes X42B10T00.
The solution to the current problem may be derived by using the solution of the analogous X24B01T00 problem obtained by Carslaw and Jaeger ([
6], p. 128, Equation (5)). In fact, it would be sufficient to replace the variable
with
. However, as the solution method is not given in [
6], a complete derivation using Green’s functions is proposed afterwards.
The Green’s function for the X42 case is available in ([
4], p. 609, Equation (42)). In dimensionless form, it results in:
where the eigenfunction,
, and the dimensionless norm,
are defined as
Additionally, the eigenvalues
satisfy the following eigencondition
whose solution is discussed in
Section 4.1.
Applying the Green’s Function Solution Equation (GFSE) ([
4], Chap. 6, pp. 181–182) to the current problem yields
Then, by substituting Equation (19a) into the above integral, it is found
whose integration gives
As the first summation exhibits an algebraic convergence (very slow), which requires a large number of terms, an alternative solution with better convergence properties is sought. This solution is discussed in
Appendix A. Thus, by using the algebraic identity defined by Equation (A15), Equation (22) becomes
Also, by using the interface condition defined by Equation (18c), the thin layer temperature results in
Equation (23) show two parts: a quasi-steady part,
or
, which yields an increasing temperature response with time, and a complementary transient part,
or
, which becomes negligible at large times. As regards the latter, an infinite number of terms cannot be taken into account. For this reason, a finite number of terms will be considered and the related truncation error (termed ‘tail’) is discussed in
Section 4.2.
4.1. Computation of the Eigenvalues
The eigenvalues appearing in the temperature solution are computed as roots of the eigencondition defined by Equation (19d). This transcendental equation is similar to the corresponding equations of the X31 = X13 cases treated in [
14]. In fact, in order to obtain Equation (19d), it is sufficient to replace the Biot number with the term
1/
P. Thus, its roots (eigenvalues) may be computed by using the same explicit approximate relations based on the second-order modified Newton method [
14]. These relations provide an approximate value of the exact eigenvalue with high accuracy (8-decimal place after one iteration, and 15-decimal place after two iterations) for the
P range
. A computer code for computing the eigenvalues of the X42 case is given in [
15].
4.2. Maximum Number of Terms
The infinite summation appearing in Equation (23) exhibits an exponential convergence and, hence, very fast. A conservative convergence criterion for this summation may be defined following the procedure given in [
16]. In particular, by using
as a conservative estimate for
, the maximum number of required terms
in order to get a truncation error of
(
A = 2, 3, …, 15) in the above series solution may be taken as:
where the function ‘ceil(
z)’ is a Matlab function which rounds the number
z to the nearest integer greater than or equal to
z.
However, the criterion defined by Equation (24) may require a large number of terms for early times. To avoid it, for short times, less than the so-called deviation time
defined as [
11,
17]
the temperature solution Equation (23) can be replaced by the semi-infinite solution of the X40B1T00 problem, Equation (3), with an error less than
.
Moreover, for large times, greater than the quasi-steady time
defined by [
15]
the complementary transient solution of Equation (23) becomes negligible with errors less than
(
A = 2, 3, …, 15). Thus, for
the temperature solution can be taken as
A verified computer code in Matlab ambient for computing the solution to the X42B10T00 problem is provided in [
15]. For the sake of completeness, some numerical values of
and
for different accuracies
A are shown in
Table 1.
4.3. Thin Layer between Two Slabs
The solution defined by Equation (23) can also be used for a three-layer configuration in which a thin layer is located between two slabs insulated at the other boundaries. This situation is depicted in
Figure 4, where the thin layer and the slabs are in perfect thermal contact and initially at uniform temperature
Tin. At
t = 0, an internal heat generation per unit of volume
occurs only within the thin layer. For the addressed problem the governing equations are expressed as
where Equation (28h,i) denote a perfect contact, and the dimensionless variables are defined by Equation (5).
For the sake of thermal symmetry, the three layer structure depicted in
Figure 4 reduces to the two-layer configuration shown in Figure 5 for perfect contact. Therefore, the thermal field
within the right slab of the three-layer configuration is still defined by Equation (23), which may also be used to describe the thermal field
of the left slab, by replacing the variable
with
. As regards the thin layer, its temperature is given by Equation (23) for
.
It is interesting to observe that the schematic of
Figure 4 represents the experimental apparatus used for thermal properties measurements [
2,
18]. In particular, the thin film represents the heater, while the two slabs are the specimens of the material under investigation.
5. Thin Layer in Imperfect Contact with a Slab
This problem is similar to the one discussed in
Section 4 and depicted in
Figure 3, but now an imperfect thermal contact with the surface layer is taken into account. Additionally, the slab and the thin layer are again at the same uniform initial temperature
Tin. The problem notation is X62B10T11.
In dimensionless form, the governing equations for the current problem (whose notation becomes X62B10T00) are still defined by Equation (18), except the interface condition, Equation (18c), which has to be replaced with Equation (6).
The exact analytical solution to the current X62B10T00 problem is available in the heat conduction literature ([
6], p. 129, Equation (14)). However, as in this reference the solution method is not given, a complete derivation using Laplace transform is proposed afterwards. Contrary to the X42 case treated in
Section 4.2, the GF approach cannot here be used as the
Green’s function is not available in the literature.
By applying the Laplace transform to Equations (18a,b,d) and (6), subject to the initial conditions, Equation (18e,f), it results in the same problem defined by Equation (7), where now Equation (7c) valid for a semi-infinite body has to be replaced with an insulating condition, that is,
where
.
In the current case however, it is convenient to rewrite the solution given by Equation (8) by using a linear combinations of exponentials, such as
and
. Thus, Equation (8) becomes
where
and the phase
can be determined by applying the boundary condition at
, i.e., Equation (29). It follows that:
. Then, as the hyperbolic cosine is an even function, Equation (30) becomes
Also, by using the boundary condition of the sixth kind defined by Equation (7b,d) in the Laplace domain, it is possible to determine the constant
A as
Thus, the sought solution
is
To take the inverse of Equation (33) a Taylor series expansion of the hyperbolic functions is needed. It results in
Therefore, the function
at the denominator of Equation (33) becomes
The above function admits a quadruple pole at
q =
0 (which implies a double pole at
s =
0). Additionally, it has simple zeros along the negative real axis, that may conveniently be assumed as
(or
). Then, the value of
can be computed by solving numerically the eigencondition
, that is
where the relations
and
have been used. Equation (36) admits an infinite number of simple poles at
,
,
m = 1, 2, … (or
), where
is the
m-th eigenvalue whose computation is discussed in
Section 5.1.
Since the function
in Equation (33) is analytic except at poles
q = 0 (
s = 0) and at poles
(
),
,
m = 1, 2, …, the inverse of such a function may be obtained using the residue theorem ([
19], pp. 384–385).
where the calculation of the residues,
and
, is given in
Appendix B. In particular, they result in
Substituting Equation (38a,b) in Equation (37) yields the temperature solution in the time domain.
where
is the dimensionless norm given by
Once
is known, the thin layer temperature
may be obtained by using the boundary condition defined by Equation (6). After some algebra, it is given by
Similarly to the X42B10T00 case, the temperature solutions given by Equation (39a,c) show two components: a quasi-steady part and a complementary transient part that becomes negligible at large times. The maximum number of terms and the truncation error of the sums in Equation (39a,c) will be analyzed in
Section 5.2.
5.1. Computation of the Eigenvalues
The X62 eigencondition defined by Equation (36) may be rewritten as
This transcendental equation is similar to the corresponding equation of the X33 case approached in [
14]. In detail, it is needed to replace the sum of the Biot numbers
with the term
, and their product
with the term
. Therefore, its roots (eigenvalues) may be computed by using the same explicit approximate relations based on the third-order modified Newton method [
14]. These relations provide an approximate value
of the exact eigenvalue
(m = 1, 2, 3, …), with at least seven-decimal place accuracy (10
−7) for
and for
. Additionally, the lowest accuracy of 10
−7 occurs when
m = 1; on the contrary, when m > 1, the accuracy is higher [
14]. A computer code for computing the eigenvalues of the X62 case is given in Ref. [
20].
5.2. Computational Solution
The exact analytical solution given by Equation (39) shows an infinite summation, while the computational analytical solution requires a finite number of terms. In particular, the maximum number of required terms
for a truncation error less than
(
A = 2, 3, …, 15), may be taken in a conservative way as [
20]:
where
denotes the Heaviside step function.
As the convergence criterion defined by Equation (41) may require a large number of terms for early times, the solution of the X62B10T00 problem can be replaced by the semi-infinite transient X60B1T00 solution given by Equations (15) and (16). In fact, at early times the finite body behaves as a semi-infinite one. In particular, for times less than the so-called deviation time
defined by Equation (25), this replacement occurs with errors with less than
[
11,
17].
Moreover, for time greater than the so-called steady time
, the complementary transient part of the solution defined by Equation (39a,c) may be neglected with errors less than
. Therefore,
where
may be taken in a conservative way as [
20]
Some numerical values of
for different accuracies
A are shown in
Table 2. Additionally, a verified computer code in Matlab ambient for computing the solution of the X62B10T00 problem is provided in [
20].
5.3. Thin Layer between Two Slabs: Imperfect Contact Case
The solution given by Equation (39) may also be used for the three-layer configuration shown in
Figure 4 when an imperfect thermal contact between the slabs and the thin layer is taken into account.
The related governing equations may still be defined by Equation (28), provided that Equation (28h,i) be replaced by the following interface conditions, respectively:
With the same argumentations reported in
Section 4.3 the addressed problem reduces to the X62B10T00 case. Hence, the thermal field within the right slab, in the three-layer configuration, is still defined by Equation (41a). It may also be used to describe the thermal field
, developing in the left slab, by replacing the variable
with
. For what concerns the thin layer, its temperature
is still given by Equation (39c).
It has already been mentioned in
Section 4.3 that the schematic of
Figure 4 finds application in the field of thermal properties measurements when using transient inverse techniques. In such an application, the contact resistance
accounts for the imperfect contact between the heater (thin layer) and the specimens (slabs) [
21].
7. Conclusions
One-dimensional multi-layer transient heat conduction problems involving boundary conditions of fourth and sixth kind have been discussed.
In particular, structures in which a thin layer and a semi-infinite body are in perfect or imperfect thermal contact (X40B1T00 and X60B1T00 problems, respectively) have been investigated. Additionally in these situations a jump in the heat flux was applied to the thin layer. Furthermore, the temperature solution of the X60B1T00 problem, which is not available in the heat conduction literature, has been derived by means of the LT method.
Two-layer problems defined by a thin layer in perfect or imperfect thermal contact with a slab, insulated at the other boundary, have also been investigated. In particular, the exact analytical solutions for the X42B10T00 and X62B10T00 problems have been derived by means of the GFSE and the LT methods, respectively.
Moreover situations in which a thin layer, with internal heat generation, is located between two adjacent domains (semi-infinite bodies or slabs) have been discussed. For sake of thermal symmetry occurring in one dimensional conductive problems, such three-layer structures have been reduced to the two-layer configurations previously analyzed.
The results have been shown that for all the considered cases the temperature levels reached in the domain increase when the heat capacity of the thin layer and the contact resistance at the thin layer-domain interface decrease. In particular, the greatest temperature field has been obtained for a finite body in perfect thermal contact with a thin layer. Additionally, when an imperfect contact is considered, it has been found that the thin layer temperature is highest when the adjacent domain is a slab and it increases when the contact resistance rises.