Next Article in Journal
Grid Code-Dependent Frequency Control Optimization in Multi-Terminal DC Networks
Next Article in Special Issue
Management of the Torch Structure with the New Methodological Approaches to Regulation Based on Neural Network Algorithms
Previous Article in Journal
Distributed Control of Clustered Populations of Thermostatic Loads in Multi-Area Systems: A Mean Field Game Approach
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Multi-Layer Transient Heat Conduction Involving Perfectly-Conducting Solids

by
Giampaolo D’Alessandro
* and
Filippo de Monte
Department of Industrial and Information Engineering and Economics, University of L’Aquila, via G. Gronchi 18, 67100 L’Aquila, Italy
*
Author to whom correspondence should be addressed.
Energies 2020, 13(24), 6484; https://doi.org/10.3390/en13246484
Submission received: 30 October 2020 / Revised: 2 December 2020 / Accepted: 5 December 2020 / Published: 8 December 2020
(This article belongs to the Special Issue Heat Diffusion: Dynamical Modelling, Control)

Abstract

:
Boundary conditions of high kinds (fourth and sixth kind) as defined by Carslaw and Jaeger are used in this work to model the thermal behavior of perfect conductors when involved in multi-layer transient heat conduction problems. In detail, two- and three-layer configurations are analyzed. In the former, a thin layer modeled as a lumped body is subject to a surface heat flux on the front side while it is in perfect (fourth kind) or in imperfect (sixth kind) thermal contact with a semi-infinite or finite body on the back side. When dealing with a semi-infinite body in imperfect contact, the temperature solution is derived by means of the Laplace transform method. Green’s function approach is also used but for solving the companion case of a finite body in perfect contact with the thin film. In the latter, a thin layer with internal heat generation is located between two semi-infinite or finite bodies in perfect/imperfect contact. For the sake of thermal symmetry, such a three-layer structure reduces to a two-layer configuration. Results are given in both tabular and graphical forms and show the effect of heat capacity and thermal resistance on the temperature distribution of conductive layers.

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 t = 0 a step change in heat flux is applied to the surface layer at its boundary x = L f , 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
2 T ˜ x ˜ 2 = T ˜ t ˜         ( 0 < x ˜ < ; t ˜ > 0 ) ,
T ˜ x ˜ | x ˜ = 0 + P T ˜ f t ˜ = 1         ( t ˜ > 0 ) ,
T ˜ f ( t ) = T ˜ ( 0 , t )         ( t ˜ > 0 ) ,
T ˜ ( x ˜ , t ˜ ) = finite           ( t ˜ > 0 ) ,
T ˜ ( x ˜ , 0 ) = 0           ( 0 < x ˜ < ) ,
T ˜ f ( 0 ) = 0
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
T ˜ = T T i n q f , 0 L / k ,   T ˜ f = T f T i n q f , 0 L / k ,   x ˜ = x L ,   t ˜ = α t L 2 ,   P = C f L f C L ,
where L is a reference length, while L f 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
T ˜ ( x ˜ , t ˜ ) = 4 t ˜ π e x ˜ 2 4 t ˜ x ˜ erfc ( x ˜ 4 t ˜ ) P [ erfc ( x ˜ 4 t ˜ ) Am ( x ˜ , t ˜ , P ) ]   ( 0 x ˜ < ; t ˜ 0 ) ,
where Am ( x ˜ , t ˜ , P ) is the so-called Amos function, defined as
Am ( x ˜ , t ˜ , P ) = e x ˜ P + t ˜ P 2 erfc ( x ˜ 4 t ˜ + t ˜ P ) = e x ˜ 2 4 t ˜ erfcx ( x ˜ 4 t ˜ + t ˜ P ) ,
The erfcx   ( z ) = e z 2 erfc   ( z ) is the scaled complementary error function. It is recommended for large values of x ˜ , t ˜ and 1 / P 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 x ˜ , t ˜ and 1 / P .

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 g f , 0 .
In such a case, the governing equations in dimensionless form are
2 T ˜ S I ( ) x ˜ 2 = T ˜ S I ( ) t ˜             ( < x ˜ < 2 L ˜ f ;   t ˜ > 0 ) ,
T ˜ S I ( ) ( x ˜ , t ˜ ) = finite               ( t ˜ > 0 ) ,
T ˜ S I ( ) ( x ˜ , 0 ) = 0             ( < x ˜ < 2 L ˜ f ) ,
2 T S I ( + ) x ˜ 2 = T S I ( + ) t ˜             ( 0 < x ˜ < ;   t ˜ > 0 ) ,
T ˜ S I ( + ) ( x ˜ , t ˜ ) = finite             ( t ˜ > 0 ) ,
T ˜ S I ( + ) ( x ˜ , 0 ) = 0                         ( 0 < x ˜ < ) ,
2 P T ˜ f t ˜ = 2 + T ˜ S I ( + ) x ˜ | x ˜ = 0 T ˜ S I ( ) x ˜ | x ˜ = 2 L ˜ f       ( t ˜ > 0 ) ,
T ˜ S I ( ) ( 2 L ˜ f , t ˜ ) = T ˜ f ( t ˜ )       ( t ˜ > 0 ) ,
T ˜ f ( t ˜ ) = T ˜ S I ( + ) ( 0 , t ˜ ) ( t ˜ > 0 ) ,
T ˜ f ( 0 ) = 0 ,
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
T ˜ = T T i n ( g f , 0     L f ) L / k ,   T ˜ f = T f T i n ( g f , 0     L f ) L / k ,   x ˜ = x L ,   t ˜ = α t L 2 ,   P = C f L f C L ,   L ˜ f = L f L ,
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 L f ), and by noting that ( g f , 0     L f ) in Equation (5) plays the same role of the surface heat flux q f , 0 in Equation (2), the problem defined by Equation (4) reduces to the X40B1T00 case defined by Equation (1). Hence, the thermal field T ˜ S I ( + ) ( x ˜ , t ˜ ) 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 T ˜ S I ( ) ( x ˜ , t ˜ ) on the left-side by simply replacing the variable x ˜ with ( x ˜ + 2 L ˜ f ) . As far as the thin film is concerned, its temperature is defined by Equation (3a) for x ˜ = 0 .

3. Thin Layer in Imperfect Contact with a Semi-Infinite Body

Consider the same semi-infinite body of the problem discussed in Section 2 and shown in Figure 1, where now an imperfect thermal contact with the surface layer at the boundary x = 0 is taken into account through a contact resistance Rc. The problem notation is X60B1T11.
The governing equations, in dimensionless form, for this conductive problem are still defined by Equation (1), except the condition in Equation (1c) which has to be replaced with the following:
T ˜ x ˜ | x ˜ = 0 = 1 R ˜ c [ T ˜ f ( t ˜ ) T ˜ ( 0 , t ˜ ) ]                         ( t ˜ > 0 ) ,
where R ˜ c = R c L / k is the dimensionless thermal contact resistance. Additionally, Equations (1b) and (6) represent the boundary condition of the sixth kind. In dimensionless form, the problem notation becomes X60B1T00.

3.1. Temperature Solution

To authors’ knowledge, the solution to the current problem is not available in the heat conduction literature. However, it may be obtained by applying Laplace Transform to Equations (1a,b,d) and (6), bearing in mind the initial conditions Equation (1e,f). It results in
2 ϑ x ˜ 2 = s ϑ                         ( 0 < x ˜ < ) ,
ϑ x ˜ | x ˜ = 0 + P s ϑ f = 1 s ,  
ϑ ( x ˜ , s ) = finite ,  
ϑ f = ϑ ( 0 , s ) R ˜ c ϑ x ˜ | x ˜ = 0 ,  
where s denotes the Laplace variable, ϑ ( x ˜ , s ) = L [ T ( x ˜ , t ˜ ) ] , and ϑ f ( s ) = L [ T ˜ f ( t ˜ ) ] .
Thus, seeking the canonical solution for Equation (7a) yields
ϑ ( x ˜ , q ) = c 1 exp   ( q   x ˜ ) + c 2 exp   ( q   x ˜ ) ,
where q = s .
The boundary condition Equation (7c) applied to Equation (8) implies c 1 = 0 . Moreover, substituting Equation (7d) into Equation (7b) gives
ϑ x ˜ | x ˜ = 0 ( 1 + P R ˜ c s ) + P s ϑ ( 0 , s ) = 1 s ,
where ϑ ( 0 , s ) and its derivative ( ϑ /   x ˜ ) x ˜ = 0 may be calculated using Equation (8). Behaving like this, the constant c 2 can be determined and the solution ϑ ( x ˜ , s ) is
ϑ ( x ˜ , q ) = exp ( q   x ˜ ) P R ˜ c q 3 ( 1 P R ˜ c + q R ˜ c + q 2 ) ,
To take the inverse of Equation (10) by means of the standard inversion tables, it is convenient to rewrite it as
ϑ ( x ˜ , q ) = exp ( q   x ˜ ) P R ˜ c     q 3     ( q + β ˜   1 )   ( q + β ˜   2 ) ,
where β ˜   1 and β ˜   2 are the roots of the polynomial appearing between round brackets at the denominator of Equation (10). They are
β ˜   1 = 1 2 R ˜ c ( 1 + 1 4 R ˜ c P ) ;   β ˜   2 = 1 2 R ˜ c ( 1 1 4 R ˜ c P ) ,
By some algebra, Equation (11) becomes
ϑ ( x ˜ , q ) = exp ( q   x ˜ ) P R ˜ c   s   q     ( q + β ˜   1 )   ( q + β ˜   2 )   ( q + β ˜   2 )   ( q + β ˜   1 )     ( β ˜   2 β ˜   1 )   = β ˜   1   β ˜   2     ( β ˜   2 β ˜   1 ) [ exp ( q   x ˜ ) s   q     ( q + β ˜   1 ) exp ( q   x ˜ ) s   q     ( q + β ˜   2 ) ]
where the relation β ˜   1   β ˜   2 = 1 / ( P R ˜ c ) has been used.
Now, taking the inverse Laplace transform of Equation (13) yields
T ˜ ( x ˜ , t ˜ ) = L   1 [ ϑ ( x ˜ , q )   ] = β ˜   1   β ˜   2     ( β ˜   2 β ˜   1 ) { L   1 [ exp ( q   x ˜ )   q 3     ( q + β ˜   1 ) ] L   1 [ exp ( q   x ˜ ) q 3     ( q + β ˜   2 ) ]   } ,
where, by using the inversion tables provided in ([6], p. 495, Equation (15)),
L   1 [ exp ( q   x ˜ ) s   q     ( q + β ˜   i ) ] = 1 β ˜   i 4 t ˜ π exp     [   x ˜ 2 / ( 4 t ˜ ) ] 1 + β ˜   i   x ˜ β ˜ i 2 erfc ( x ˜ 4 t ˜ ) + 1 β ˜ i 2 Am ( x ˜ , t ˜ , 1 / β ˜   i )             ( i   =   1 , 2 )
Note that Am ( x ˜ , t ˜ , 1 / β ˜   i ) is the Amos function defined by Equation (3b) for P = 1 / β ˜   i .
Then, by substituting Equation (14b) into Equation (14a), the temperature solution is
T ˜ ( x ˜ , t ˜ ) = 4 t ˜ π exp     [   x ˜ 2 / ( 4 t ˜ ) ] ( P + x ˜ ) erfc ( x ˜ 4 t ˜ ) + β ˜   2   ( β ˜   2 β ˜   1 )   β ˜   1 Am ( x ˜ , t ˜ , 1 / β ˜   1 ) β ˜   1   ( β ˜   2 β ˜   1 )   β ˜   2 Am ( x ˜ , t ˜ , 1 / β ˜   2 )             ( 0 x ˜ < ;   t ˜ 0 ) ,
Once the semi-infinite solution is known, the surface layer temperature may be determined through the boundary condition defined by Equation (6). It is
T ˜ f ( t ˜ ) = 4 t ˜ π P + R ˜ c + 1   ( β ˜   2 2 β ˜   1 2 )   [ β ˜   2 2 β ˜   1 Am ( 0 , t ˜ , 1 / β ˜   1 ) β ˜   1 2 β ˜   2 Am ( 0 , t ˜ , 1 / β ˜   2 ) ]       ( t ˜ 0 ) ,
The solution given by Equation (15) is valid for any value of ( 4 R ˜ c / P ) greater than, equal or less than 1. Nevertheless, when 4 R ˜ c / P = 1 , it is convenient to define the solution as follows
T ˜ ( x ˜ , t ˜ ) = 2 4 t ˜ π e x ˜ 2 4 t ˜ ( P + x ˜ ) erfc ( x ˜ 4 t ˜ ) + ( P x ˜ 4 t ˜ P ) Am ( x ˜ , t ˜ , P / 2 ) ,
T ˜ f ( t ˜ ) = 3 2 4 t ˜ π 3 4 P + ( 3 4 P 2 t ˜ P )   erfcx ( 2 t ˜ P ) ,
A verified computer code for calculating the solution of the X60B1T00 problem is provided in [13].

3.2. Thin Layer between Two Semi-Infinite Bodies: Imperfect Contact Case

The solution given by Equations (15) and (16) can also be used for the three-layer configuration shown in Figure 2 when an imperfect thermal contact between the thin layer and the semi-infinite bodies occurs. In such a case, the governing equations may still be defined by Equation (4), provided that Equation (4h,i) be replaced by the following conditions, respectively:
T ˜ S I ( ) x ˜ | x ˜ = 2 L ˜ f = 1 R ˜ c [ T ˜ S I ( ) ( 2 L ˜ f , t ˜ ) T ˜ f ( t ˜ ) ]             ( t ˜ > 0 ) ,
T ˜ S I ( + ) x ˜ | x ˜ = 0 = 1 R ˜ c [ T ˜ f ( t ˜ ) T ˜ S I ( + ) ( 0 , t ˜ ) ]             ( t ˜ > 0 ) ,
where the dimensionless variables are still defined by Equation (5), and R ˜ c = R c / ( L / k ) .
In this case, by virtue of the same considerations reported in Section 2.1, the addressed problem reduces to the X60B1T00 case. Hence, the thermal field T ˜ S I ( + ) ( x ˜ , t ˜ ) within the semi-infinite body on the right-side of the three-layer configuration is still defined by Equations (15a) and (16a). For what concerns the thermal field T ˜ S I ( ) ( x ˜ , t ˜ ) , on the left-side of Figure 2, Equations (15a) and (16a) may still be used by replacing the variable x ˜ with ( x ˜ + 2 L ˜ f ) . As regards the surface layer, its temperature T ˜ f ( t ˜ ) is given by Equations (15b) and (16b).

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 t = 0 a step change in heat flux is applied to the thin film at x = L f , 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.
2 T ˜ x ˜ 2 = T ˜ t ˜       ( 0 < x ˜ < 1 ; t ˜ > 0 ) ,  
T ˜ x ˜ | x ˜ = 0 + P T ˜ f t ˜ = 1       ( t ˜ > 0 ) ,
T ˜ f ( t ) = T ˜ ( 0 , t )       ( t ˜ > 0 ) ,  
( T ˜ x ˜ ) x ˜ = 1 = 0             ( t ˜ > 0 ) ,
T ˜ ( x ˜ , 0 ) = 0     ( 0 < x ˜ < 1 ) ,
T ˜ f ( 0 ) = 0 ,  
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 x ˜ with 1 x ˜ . 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:
G ˜ X 42 ( x ˜ , t ˜ | x ˜ , τ ˜ ) = 1 P + 1 + 2 m = 1 X m ( x ˜ ) X m ( x ˜ ) N ˜ m exp [ β m 2 ( t ˜ τ ˜ ) ] ,
where the eigenfunction, X m ( x ˜ ) , and the dimensionless norm, N ˜ m are defined as
X m ( x ˜ ) = cos ( β m x ˜ ) P β m sin ( β m x ˜ ) ,
N ˜ m = ( P β m ) 2 + P + 1 ,
Additionally, the eigenvalues β m satisfy the following eigencondition
β m cot ( β m ) = 1 P ,
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
T ˜ ( x ˜ , t ˜ ) = τ ˜ = 0 t ˜ G ˜ X 42 ( x ˜ , t ˜ | 0 , τ ˜ ) d τ ˜ ,
Then, by substituting Equation (19a) into the above integral, it is found
T ˜ ( x ˜ , t ˜ ) = τ ˜ = 0 t ˜ 1 P + 1 d τ ˜ + 2 m = 1 X m ( x ˜ ) N ˜ m τ ˜ = 0 t ˜ e β m 2 ( t ˜ τ ˜ ) d τ ˜ ,
whose integration gives
T ˜ ( x ˜ , t ˜ ) = t ˜ P + 1 + 2 m = 1 cos ( β m x ˜ ) P β m sin ( β m x ˜ ) N ˜ m β m 2 2 m = 1 cos ( β m x ˜ ) P β m ( β m x ˜ ) N ˜ m β m 2 exp ( β m 2 t ˜ ) ,
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
T ˜ ( x ˜ , t ˜ ) = t ˜ P + 1 + x ˜ 2 2 ( P + 1 ) x ˜ P + 1 + 1 3 ( P + 1 ) 2 T ˜ s     ( x ˜ , t ˜ ) 2 m = 1 cos ( β m x ˜ ) P β m sin ( β m x ˜ ) N ˜ m β m 2 e β m 2 t ˜ T ˜ c t     ( x ˜ , t ˜ ) ( 0 x ˜ 1 ; t ˜ 0 ) ,
Also, by using the interface condition defined by Equation (18c), the thin layer temperature results in
T ˜ f   ( t ˜ ) = t ˜ P + 1 + 1 3 ( P + 1 ) 2 T ˜ f , s     ( t ˜ ) 2 m = 1 e β m 2 t ˜ N ˜ m β m 2 T ˜ f , c t     ( t ˜ )             ( t ˜ 0 ) ,
Equation (23) show two parts: a quasi-steady part, T ˜ s ( x ˜ , t ˜ ) or T ˜ f ,   s ( t ˜ ) , which yields an increasing temperature response with time, and a complementary transient part, T ˜ c t ( x ˜ , t ˜ ) or T ˜ f ,   c t ( t ˜ ) , 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 [ 0 , ) . 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 ( m 1 / 2 ) π as a conservative estimate for β m , the maximum number of required terms m max in order to get a truncation error of 10 A (A = 2, 3, …, 15) in the above series solution may be taken as:
m max = c e i l { 1 2 + 1 π [ A ln ( 10 ) t ˜ ] 1 / 2 } ,
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 t ˜ d defined as [11,17]
t ˜ d = ( 2 x ˜ ) 2 / ( 10 A )                         ( A =   2 ,   3 ,   ,   15 ) ,
the temperature solution Equation (23) can be replaced by the semi-infinite solution of the X40B1T00 problem, Equation (3), with an error less than 10 A .
Moreover, for large times, greater than the quasi-steady time t ˜ s defined by [15]
t ˜ s = 4 π 2 { A ln 10 + ln [ 24 ( P + 1 ) 2 π 2 ] } ,
the complementary transient solution of Equation (23) becomes negligible with errors less than 10 A (A = 2, 3, …, 15). Thus, for t ˜ t ˜ s the temperature solution can be taken as
T ˜ ( x ˜ , t ˜ ) = t ˜ P + 1 + x ˜ 2 2 ( P + 1 ) x ˜ P + 1 + 1 3 ( P + 1 ) 2                   ( 0 x ˜ 1 ; t ˜ t ˜ s ) ,
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 t ˜ d and t ˜ s 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 g f , 0 occurs only within the thin layer. For the addressed problem the governing equations are expressed as
2 T ˜ s l a b ( ) x ˜ 2 = T ˜ s l a b ( ) t ˜             ( ( 1 + 2 L ˜ f ) < x ˜ < 2 L ˜ f ; t ˜ > 0 ) ,
T ˜ s l a b ( ) x ˜ | x ˜ = ( 1 + 2 L ˜ f ) = 0             ( t ˜ > 0 ) ,
T ˜ s l a b ( ) ( x ˜ , 0 ) = 0           ( ( 1 + 2 L ˜ f ) < x ˜ < 2 L ˜ f ) ,
2 T s l a b ( + ) x ˜ 2 = T s l a b ( + ) t ˜             ( 0 < x ˜ < 1 ; t ˜ > 0 ) ,
T ˜ s l a b ( + ) x ˜ | x ˜ = 1 = 0             ( t ˜ > 0 ) ,
T ˜ s l a b ( + ) ( x ˜ , 0 ) = 0           ( 0 < x ˜ < 1 ) ,
2 P T ˜ f t ˜ = 2 + T ˜ s l a b ( + ) x ˜ | x ˜ = 0 T ˜ s l a b ( ) x ˜ | x ˜ = 2 L ˜ f             ( t ˜ > 0 ) ,
T ˜ s l a b ( ) ( 2 L ˜ f , t ˜ ) = T ˜ f ( t ˜ )             ( t ˜ > 0 ) ,
T ˜ f ( t ˜ ) = T ˜ s l a b ( + ) ( 0 , t ˜ )             ( t ˜ > 0 ) ,
T ˜ f ( 0 ) = 0 ,
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 T ˜ s l a b ( + ) ( x ˜ , t ˜ ) 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 T ˜ s l a b ( ) ( x ˜ , t ˜ ) of the left slab, by replacing the variable x ˜ with ( x ˜ + 2 L ˜ f ) . As regards the thin layer, its temperature is given by Equation (23) for x ˜ = 0 .
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 G X 62 ( x ˜ , t ˜ | x ˜ , τ ˜ ) 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,
( ϑ x ˜ ) x ˜ = 1 = 0 ,
where ϑ ( x ˜ , s ) = L [ T ˜ ( x ˜ , t ˜ ) ] .
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 ( e q x ˜ e q x ˜ ) / 2 and ( e q x ˜ + e q x ˜ ) / 2 . Thus, Equation (8) becomes
ϑ ( x ˜ , q ) = B   1 cosh ( q x ˜ ) + B   2 sinh ( q x ˜ ) = A cosh ( q x ˜ + ϕ ) ,
where q = s and the phase ϕ can be determined by applying the boundary condition at x ˜ = 1 , i.e., Equation (29). It follows that: ϕ = q . Then, as the hyperbolic cosine is an even function, Equation (30) becomes
ϑ ( x ˜ , q ) = A cosh [ q ( x ˜ 1 ) ] = A cosh [ q ( 1 x ˜ ) ] ,
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
A = 1 P q 4 cosh ( q ) + q   3 ( P R ˜ c   q 2 + 1 ) sinh ( q ) ,
Thus, the sought solution ϑ ( x ˜ , q ) is
ϑ ( x ˜ , q ) = cosh [ q ( 1 x ˜ ) ] q 2 [ P q 2 cosh ( q ) + q ( P R ˜ c   q 2 + 1 ) sinh ( q ) ] D ( q ) = cosh [ q ( 1 x ˜ ) ] D ( q ) ,
To take the inverse of Equation (33) a Taylor series expansion of the hyperbolic functions is needed. It results in
cosh [ q ( 1 x ˜ ) ] = 1 + q 2 ( 1 x ˜ ) 2 + O ( q 4 ) ,
cosh ( q ) = 1 + q 2 2 + O ( q 4 ) ,
sinh ( q ) = q + q 3 3 ! + O ( q 5 ) ,
Therefore, the function D ( q ) at the denominator of Equation (33) becomes
D ( q ) = q 2 { P q 2 [ 1 + q 2 2 + O ( q 4 ) ] + q ( P R ˜ c   q 2 + 1 ) [ q + q 3 3 ! + O ( q 5 ) ] } = q 4 { P + 1 + q 2 [ P 2 + P R ˜ c + 1 6 ] + O ( q 4 ) }
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 q = i β (or s = β 2 ). Then, the value of β can be computed by solving numerically the eigencondition D ( q = i β ) = 0 , that is
P β cos ( β ) + ( 1 P R ˜ c   β 2 ) sin ( β ) = 0 ,
where the relations cosh ( i β ) = cos ( β ) and sinh ( i β ) = i sin ( β ) have been used. Equation (36) admits an infinite number of simple poles at q = i β m ,   β m > 0 , m = 1, 2, … (or s = β m 2 ), where   β m is the m-th eigenvalue whose computation is discussed in Section 5.1.
Since the function ϑ ( x ˜ , q ) in Equation (33) is analytic except at poles q = 0 (s = 0) and at poles q = i β m ( s = β m 2   ), β m > 0 , m = 1, 2, …, the inverse of such a function may be obtained using the residue theorem ([19], pp. 384–385).
T ˜ ( x ˜ , t ˜ ) = L 1 { ϑ ( x ˜ , s ) } = 1 2 π   i γ i γ + i ϑ ( x ˜ , s )     e s t d s =   Residues in the closed contour     encompassing the left half plane = R s = 0 + m = 1 R s = β m 2
where the calculation of the residues, R s = 0 and R s = β m 2 , is given in Appendix B. In particular, they result in
R s = 0 = t ˜ P + 1 + x ˜ 2 2 ( P + 1 ) x ˜ P + 1 + 1 3 P R ˜ c 3 ( P + 1 ) 2 ,
R s = β m 2 = 2   ( P R ˜ c β m 2 1 ) cos [ β m ( 1 x ˜ ) ] e β m 2 t ˜ [ ( P R ˜ c ) 2 β m 4 + ( P + P R ˜ c 2 R ˜ c ) P β m 2 + P + 1 ] β m 2 cos ( β m ) ,
Substituting Equation (38a,b) in Equation (37) yields the temperature solution in the time domain.
T ˜ ( x ˜ , t ˜ ) = t ˜ P + 1 + x ˜ 2 2 ( P + 1 ) x ˜ P + 1 + 1 3 P R ˜ c 3 ( P + 1 ) 2 T ˜ s ( x ˜ , t ˜ ) + 2 m = 1 ( P R ˜ c β m 2 1 ) cos [ β m ( 1 x ˜ ) ] N ˜ m β m 2 cos ( β m ) e β m 2 t ˜ T ˜ c t ( x ˜ , t ˜ ) ( 0 x ˜ 1 ; t ˜ 0 ) ,
where N ˜ m is the dimensionless norm given by
N ˜ m = ( P R ˜ c ) 2 β m 4 + ( P + P R ˜ c 2 R ˜ c ) P β m 2 + P + 1 ,
Once T ˜ ( x ˜ , t ˜ ) is known, the thin layer temperature T ˜ f ( t ˜ ) may be obtained by using the boundary condition defined by Equation (6). After some algebra, it is given by
T ˜ f ( t ˜ ) = t ˜ P + 1 + R ˜ c + 1 / 3 ( P + 1 ) 2 T ˜ f , s (   t ˜ ) 2 m = 1 e β m 2 t ˜ N ˜ m β m 2 T ˜ f , c t (   t ˜ )                   ( t ˜ 0 ) ,
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
tan ( β m ) = β m R ˜ c β m   2 1 P R ˜ c ,
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 ( B 1 + B 2 ) with the term 1 / R ˜ c , and their product ( B 1 B 2 ) with the term 1 / ( P R ˜ c ) . 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 ζ m of the exact eigenvalue β m (m = 1, 2, 3, …), with at least seven-decimal place accuracy (10−7) for P ( 0 , ) and for R ˜ c ( 0 , ) . 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 m max for a truncation error less than 10 A (A = 2, 3, …, 15), may be taken in a conservative way as [20]:
m max = 1 + ceil { 1 π A ln ( 10 ) + [ H ( P 1 ) 1 ] ln ( P ) t ˜ } ,
where H ( ) 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 t ˜ d defined by Equation (25), this replacement occurs with errors with less than 10 A [11,17].
Moreover, for time greater than the so-called steady time t ˜ s , the complementary transient part of the solution defined by Equation (39a,c) may be neglected with errors less than 10 A . Therefore,
T ˜ ( x ˜ , t ˜ ) = t ˜ P + 1 + x ˜ 2 2 ( P + 1 ) x ˜ P + 1 + 1 3 P R ˜ c 3 ( P + 1 ) 2                   ( 0 x ˜ 1 ; t ˜ t ˜ s ) ,
T ˜ f ( t ˜ ) = t ˜ P + 1 + R ˜ c + 1 / 3 ( P + 1 ) 2                   ( t ˜ t ˜ s ) ,
where t ˜ s may be taken in a conservative way as [20]
t ˜ s = A ( 1 + 2 / A ) ln 10 + [ H ( N ˜ 1   β 1 2 1 ) 1 ] ln ( N ˜ 1   β 1 2 ) β 1 2 ,
Some numerical values of t ˜ s 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:
T ˜ s l a b ( ) x ˜ | x ˜ = 2 L ˜ f = 1 R ˜ c [ T ˜ s l a b ( ) ( 2 L ˜ f , t ˜ ) T ˜ f ( t ˜ ) ]             ( t ˜ > 0 ) ,
T ˜ s l a b ( + ) x ˜ | x ˜ = 0 = 1 R ˜ c [ T ˜ f ( t ˜ ) T ˜ s l a b ( + ) ( 0 , t ˜ ) ]                   ( t ˜ > 0 ) ,
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 T ˜ s l a b ( ) ( x ˜ , t ˜ ) , developing in the left slab, by replacing the variable x ˜ with ( x ˜ + 2 L ˜ f ) . For what concerns the thin layer, its temperature T ˜ f ( t ˜ ) 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 R ˜ c accounts for the imperfect contact between the heater (thin layer) and the specimens (slabs) [21].

6. Results and Discussion

6.1. Thin Layer in Perfect Contact with a Semi-Infinite Body

Figure 5a shows the dimensionless temperature for the X40B1T00 problem as a function of x ˜ and for different t ˜ values, while Figure 5b depicts the dimensionless temperature for the same case as a function of t ˜ and for different dimensionless locations x ˜ . Both figures are obtained for P = 1. Note that in Figure 5b the curve for x ˜ = 0 also represents the thin layer temperature T ˜ f ( t ˜ ) .
Figure 6 shows the temperature distribution within the semi-infinite medium for different values of the heat capacity ratio P and for t ˜ = 0.2 . As shown by the figure the temperature reached in the medium decreases when the P parameter rises. Additionally, it is interesting to observe that for P 0 the X40B1T00 problem reduces to the case of a semi-infinite body subject to a jump in heat flux (X20B1T0 case [22]).

6.2. Thin Layer in Imperfect Contact with a Semi-Infinite Body

Figure 7a,b depict the dimensionless temperature for the X60B1T00 problem as a function of dimensionless space and time t ˜ , respectively (for P = 1 and R ˜ c = 0.25 ). Moreover, in Figure 7b the thin layer temperature T ˜ f ( t ˜ ) is also plotted.
A comparison between Figure 5 and Figure 7 exhibits that the temperature reached for the X60B1T00 problem (imperfect thermal contact case) is slightly lower than that obtained when the perfect contact at the thin layer/semi-infinite body interface is considered (X40B1T00 problem).
Figure 8 displays the temperature distribution within the semi-infinite body for different values of the dimensionless contact resistance R ˜ c , for P = 0.1 and t ˜ = 0.2 . It can be observed that the temperature field decreases when the contact resistance increases. However, the temperature at x ˜ = 1 is less affected by the R ˜ c parameter. Note that for R ˜ c 0 the X60B1T00 problem simply reduces to the X40B1T00 problem treated in Section 2.

6.3. Thin Layer in Perfect Contact with a Slab

The dimensionless temperature field within the slab of the X42B10T00 problem (P = 1) is plotted as a function of x ˜ in Figure 9a, and as a function of time in Figure 9b. Note that in Figure 9b the curve for x ˜ = 0 also represents the thin layer temperature T ˜ f ( t ˜ ) .
By using the thickness of the slab L as a reference length for the semi-infinite body, it is possible to compare the results obtained for the X40B1T00 (Figure 5) and X42B10T00 (Figure 9) problems. Therefore, a comparison between Figure 5 and Figure 9 gives evidence that the temperature distribution within the finite slab is greater than that obtained for the semi-infinite body (not for times less than the deviation time). However, for small values of x ˜ this difference is more visible at large t ˜ values, while for large values of x ˜ it can be observed also for shorter times.
The influence of the heat capacity ratio P on the temperature distribution occurring within the slab is shown by Figure 10. A behavior similar to that displayed by Figure 6 can be observed. Note that for P 0 the X42B10T00 problem reduces to the case of a finite slab subject to a jump in heat flux (X22B10T0 case [23]).

6.4. Thin Layer in Imperfect Contact with a Slab

Figure 11 shows a plot of the dimensionless temperature as a function of x ˜ for different time values (Figure 11a), and as a function of t ˜ for different locations (Figure 11b). Additionally, the thin layer temperature T ˜ f ( t ˜ ) is depicted in Figure 11b.
By comparing Figure 9 (perfect contact case) and Figure 11 it is possible to state that the temperature field obtained for a slab in imperfect contact with a thin layer is slightly lower than that occurring when a perfect contact is considered (X42B10T00 problem).
Also, with the same argumentation reported in the previous Subsection, Figure 7 (X60B1T00 problem) and Figure 11 are compared. By observing these figures it is clear that the finite thickness of the slab yields a temperature distribution greater than that obtained in a semi-infinite body in imperfect contact with a thin layer (Figure 7). This behavior is more evident for large time values and occurs also for the thin layer temperature.
Finally, the influence of the contact resistance R ˜ c on the temperature distribution within the slab (at t ˜ = 0.2 ) is shown by Figure 12. Similarly to what has been discussed for Figure 8, when the contact resistance increases, lower temperature values are obtained. Note that for R ˜ c 0 the X62B10T00 problem simply reduces to the X42B10T00 problem discussed in Section 5.

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.

Author Contributions

Conceptualization, G.D. and F.d.M.; methodology, G.D. and F.d.M.; software, G.D.; investigation, G.D. and F.d.M.; writing—original draft preparation, G.D. and F.d.M; writing—review and editing, G.D. and F.d.M.; visualization, G.D. and F.d.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

In this Section the derivation of an alternative solution for the X42B10T00 problem is provided by means of the Alternative GFSE (AGFSE) [4].
Firstly the temperature solution is rewritten as the sum of two terms:
T ˜ ( x ˜ , t ˜ ) = T ˜ * ( x ˜ , t ˜ ) + T ˜ ( x ˜ , t ˜ ) ,
where T ˜ * ( x ˜ , t ˜ ) is the quasi-steady solution, which satisfies only the boundary conditions, and T ˜ ( x ˜ , t ˜ ) represents the transient solution.
Then, by substituting Equation (A1) into the governing equations, Equation (18), the original problem can be split up into a set of two simpler problems:
  • a non-homogeneous problem defined for T ˜ * ( x ˜ , t ˜ ) , as:
    2 T ˜ * x ˜ 2 = T ˜ * t ˜ ,
    T ˜ * x ˜ | x ˜ = 0 + P T ˜ * ( 0 , t ˜ ) t ˜ = 1 ; T ˜ * x ˜ | x ˜ = 1 = 0 ,
  • and a transient problem defined for T ˜ ( x ˜ , t ˜ ) , with homogeneous boundary conditions:
    2 T ˜ x ˜ 2 = T ˜ t ˜ ,
    T ˜ x ˜ | x ˜ = 0 + P T ˜ ( 0 , t ˜ ) t ˜ = 0 ; T ˜ x ˜ | x ˜ = 1 = 0 ,
    T ˜ ( x ˜ , 0 ) = T ˜ * ( x ˜ , 0 ) ,
    Note that in Equation (A2) the T ˜ * solution does not need to satisfy the initial condition. In fact, as the temperature solution defined by Equation (22), contains a time dependent term which dominates temperature for large times, that is t ˜ / ( P + 1 ) , the T ˜ * ( x ˜ , t ˜ ) solution should be:
    T ˜ * ( x ˜ , t ˜ ) = f ( x ˜ ) + t ˜ P + 1 ,
    where the f ( x ˜ ) function must satisfy the boundary conditions of the problem defined by Equation (A2). By substituting Equation (A4) in Equation (A2) the first problem becomes:
    d 2 f ( x ˜ ) d x ˜ 2 = 1 P + 1 ,
    d   f ( x ˜ ) d x ˜ | x ˜ = 0 + P P + 1 = 1   ,
      d   f ( x ˜ ) d x ˜ | x ˜ = 1 = 0 ,
By integrating twice Equation (A5a), the f ( x ˜ ) function is obtained:
f ( x ˜ ) = x ˜ 2 2 ( P + 1 ) + C 1 x ˜ + C 2 ,
The constants appearing in Equation (A6) can be determined by imposing the boundary conditions defined by Equation (A5b,c). In particular since the boundary conditions are both gradient conditions, the choice of the constant C2 is arbitrary. For this reason, the constant C2 is set equal to zero. Then
f ( x ˜ ) = x ˜ 2 2 ( P + 1 ) x ˜ ( P + 1 ) ,
Thus the T ˜ * ( x ˜ , t ˜ ) solution results in
T ˜ * ( x ˜ , t ˜ ) = t ˜ P + 1 + x ˜ 2 2 ( P + 1 ) x ˜ ( P + 1 ) ,
Now the transient problem defined by Equation (A3) is conveniently solved by means of the GFSE: in this case the only non-homogeneous term is the initial condition, Equation (A3c). The solution T ˜ ( x ˜ , t ˜ ) may be obtained by the following equation:
T ˜ ( x ˜ , t ˜ ) = x ˜ = 0 1 G ˜ X 42 ( x ˜ , t ˜ | x ˜ , 0 ) f ( x ˜ ) d x ˜ P G ˜ X 42 ( x ˜ , t ˜ | 0 , 0 ) f ( 0 ) = 0 ,
In the above equation the relation T ˜ * ( x ˜ , 0 ) = f ( x ˜ ) , derived from Equation (A4), has been used. Substituting Equations (19) and (A7) into the above equation gives:
T ˜ ( x ˜ , t ˜ ) = 1 P + 1 x ˜ = 0 1 ( x ˜ 2 2 ( P + 1 ) x ˜ P + 1 ) d x ˜ 2 m = 1 X m ( x ˜ ) N ˜ m e β m 2 t ˜ x ˜ = 0 1 X m ( x ˜ ) ( x ˜ 2 2 ( P + 1 ) x ˜ P + 1 ) d x ˜ ,
Bearing in mind the definition of the eigenfunction X m ( x ˜ ) , Equation (A10) becomes
T ˜ ( x ˜ , t ˜ ) = 1 2 ( P + 1 ) 2 0 1 x ˜ 2 d x ˜ + 1 ( P + 1 ) 2 0 1 x ˜ d x ˜ 2 m = 1 X m ( x ˜ ) N ˜ m e β m 2 t ˜ { 1 2 ( P + 1 ) 0 1 x ˜ 2 cos ( β m x ˜ )     d x ˜ 1 P + 1 0 1 x ˜ cos ( β m x ˜ )     d x ˜ P β m 2 ( P + 1 ) 0 1 x ˜ 2 sin ( β m x ˜ )     d x ˜ + P β m P + 1 0 1 x ˜ sin ( β m x ˜ )     d x ˜ } ,
Then by performing the integrals by parts, and after some algebra, Equation (A11) yields
T ˜ ( x ˜ , t ˜ ) = 1 3 ( P + 1 ) 2 2 m = 1 X m ( x ˜ ) N ˜ m e β m 2 t ˜ [ sin β m ( P + 1 ) β m sin β m ( P + 1 ) β m 3 P cos β m P + 1 P cos β m ( P + 1 ) β m 2 + 1 β m 2 ] ,
which, by using the eigencondition defined by Equation (19d), becomes
T ˜ ( x ˜ , t ˜ ) = 1 3 ( P + 1 ) 2 2 m = 1 X m ( x ˜ ) N ˜ m β m 2 e β m 2 t ˜ = 1 3 ( P + 1 ) 2 2 m = 1 cos ( β m x ˜ ) P β m sin ( β m x ˜ ) N ˜ m β m 2 exp ( β m 2 t ˜ ) ,
Substituting the expressions obtained for T ˜ * and T ˜ solutions, Equations (A8) and (A13), respectively, into Equation (A1), gives
T ˜ ( x ˜ , t ˜ ) = t ˜ P + 1 + x ˜ 2 2 ( P + 1 ) x ˜ P + 1 + 1 3 ( P + 1 ) 2 2 m = 1 cos ( β m x ˜ ) P β m sin ( β m x ˜ ) N ˜ m β m 2 e β m 2 t ˜ ,
It is worth noting that by setting Equations (22) and (A14) equal, the following algebraic identity can be obtained.
m = 1 cos ( β m x ˜ ) P β m sin ( β m x ˜ ) N ˜ m β m 2 = x ˜ 2 4 ( P + 1 ) x ˜ 2 ( P + 1 ) + 1 6 ( P + 1 ) 2             ( 0 x ˜ 1 ) ,
Therefore, the exact value of the infinite series has been found.

Appendix B

In this Section the calculation of the residues R s = 0 and R s = β m 2 , related to the X62B10T00 problem, is reported.
In order to determine the residue at q = 0 (or s = 0), the inversion integrand ϑ ( x ˜ , q ) e q 2   t ˜ (or ϑ ( x ˜ , s ) e s t ˜ ) should be brought into the form of a Laurent expansion at q = 0; the sought residue corresponds to the coefficient of the 1 / q 2 ( 1 / s ) term appearing in the expansion.
By using a Taylor expansion for e q 2   t ˜ , that is e q 2   t ˜ = 1 + t ˜ q 2 + t ˜ 2 2 q 4 + , and bearing in mind that ϑ ( x ˜ , q ) is a quotient of even powers of q, as shown by Equations (34) and (35), the inversion integrand can be rewritten as:
ϑ ( x ˜ , q ) e q 2   t ˜ = C 0 + C   1   q 2 + O ( q 4 ) q 4 1 + t ˜ q 2 + O ( q 4 ) C 3 + C   4   q 2 + O ( q 4 ) ,
which may be taken in a form closer to the Laurent expansion
ϑ ( x ˜ , q ) e q t ˜ = ( C 0 q 4 + C   1   q 2 + O ( 1 ) ) [ C 6 + C 7   q 2 + O ( q 4 ) ] ,
Thus, bearing in mind Equations (34a) and (35), and by setting Equation (A16a,b) equal, one can obtain the following identities:
C 0 + C   1   q 2 + O ( q 4 ) = 1 + q 2 ( 1 x ˜ ) 2 + O ( q 4 ) ,
q 4 [ C 3 + C   4   q 2 + O ( q 4 ) ] = q 4 { P + 1 + q 2 [ P 2 + P R ˜ c + 1 6 ] + O ( q 4 ) } ,
1 + t ˜ q 2 + O ( q 4 ) C 3 + C   4   q 2 + O ( q 4 ) = C 6 + C 7   q 2 + O ( q 4 ) ,
Then Equation (A17a,b) give:
C 0   = 1 , C   1 = ( 1 x ˜ ) 2 2 , C   3 = P + 1 , C   4 = P 2 + P R ˜ c + 1 6 ,
also, Equation (A17c) is an identity and the constants C   6 and C   7 can be obtained similarly.
C   6 = 1 C 3 = 1 P + 1 , C   7 = 1 C 3 ( t ˜ C 4 C 3 ) ,
Moreover, by using Equation (A16b), the coefficient of 1 / q 2 ( 1 / s ) term in the Laurent expansion may be identified. This coefficient represents the residue ( R s = 0 ) for the pole s = 0, and it results in
R s = 0 = C 0 C 7 + C 1 C 6 ,
Then by substituting Equations (A18) and (A19) into Equation (A20), and after some algebra, the same expression reported in Equation (38a) can be obtained.
Now consider the generic pole at q = i β m (or s = β m 2 ). As first step, it is worth noting that the inversion integrand is a quotient of a numerator and a denominator, i.e., ϑ ( x ˜ , q ) e q 2   t ˜ = N ( q ) / D ( q ) ; in such a case the residue for the simple pole q = i β m may be conveniently calculate by taking the derivative of the denominator of Equation (33) evaluated at its zero, that is
R s = β m 2 = N ( q = i β m ) D ( q = i β m ) = cosh [ q ( 1 x ˜ ) ]   e q 2 t ˜ | q = i β m d D ( q ) d s | q = i β m ,
where D ( q ) = q 2 [ P q 2 cosh ( q ) + q ( P R ˜ c   q 2 + 1 ) sinh ( q ) ] , see Equation (35), and its derivative d D ( q ) / d s is calculated by using the chain rule as follows
d D ( q ) d s = d D ( q ) d q d q d s = 1 2 q d D ( q ) d q ,
Thus, substituting Equation (A22) into Equation (A21) gives
R s = β m 2 = cosh [ q ( 1 x ˜ ) ] e q 2 t ˜ | q = i β m 1 2 q d D ( q ) d q | q = i β m ,
By performing the derivative appearing in Equation (A23), using the eigencondition in the form sin ( β ) = P β cos ( β ) / ( 1 P   R ˜ c   β 2 ) , after lengthy algebra, one can obtain:
1 2 q d D d q | q = i β m = 1 2 {   [ P R ˜ c β m 4 ( 4 P + 1 ) β m 2 ] cos ( β m ) + [ ( 1 + 5 R ˜ c ) P β m 3 3 β m ] sin ( β m ) } = [ ( P R ˜ c ) 2 β m 4 + ( P + P R ˜ c 2 R ˜ c ) P β m 2 + P + 1 ]     β m 2 cos ( β m ) 2 ( P R ˜ c β m 2 1 ) ,
Therefore, substituting Equation (A24) into Equation (A23) yields the residue for the generic pole s = β m 2
R s = β m 2 = 2   ( P R ˜ c β m 2 1 ) cos [ β m ( 1 x ˜ ) ] e β m 2 t ˜ [ ( P R ˜ c ) 2 β m 4 + ( P + P R ˜ c 2 R ˜ c ) P β m 2 + P + 1 ] β m 2 cos ( β m ) ,
which is exactly the same as Equation (38b).

References

  1. Van der Tempel, L.; Potze, W.; Lammers, J.H. Transient heat conduction in a thin layer between semi-infinite media in polymer shaping. J. Heat Transf. 2018, 140, 041301. [Google Scholar] [CrossRef]
  2. D’Alessandro, G.; de Monte, F. Optimal experiment design for thermal property estimation using a boundary condition of the fourth kind with a time-limited heating period. Int. J. Heat Mass Transf. 2019, 134, 1268–1282. [Google Scholar] [CrossRef]
  3. Carr, E.J.; Wood, C.J. Rear-surface integral method for calculating thermal diffusivity: Finite pulse time correction and two-layer samples. Int. J. Heat Mass Transf. 2019, 144, 118609. [Google Scholar] [CrossRef] [Green Version]
  4. Cole, K.D.; Beck, J.V.; Haji-Sheikh, A.; Litkouhi, B. Heat Conduction Using Green’s Function, 2nd ed.; CRC Press Taylor & Francis: Boca Raton, FL, USA, 2011. [Google Scholar]
  5. Al-Nimr, M.A.; Alkam, M.K. A generalized thermal boundary condition. Heat Mass Transf. 1997, 33, 157–161. [Google Scholar] [CrossRef]
  6. Carslaw, H.S.; Jaeger, J.C. Conduction of Heat in Solids, 2nd ed.; Oxford University Press: New York, NY, USA, 1959. [Google Scholar]
  7. Deconinck, B.; Pelloni, B.; Sheils, N.E. Non-steady-state heat conduction in composite walls. Proc. R. Soc. A 2014, 470, 20130605. [Google Scholar] [CrossRef] [PubMed]
  8. Rodrigo, M.R.; Worthy, A.L. Solution of multilayer diffusion problems via the Laplace transform. J. Math. Anal. Appl. 2016, 444, 475–502. [Google Scholar] [CrossRef]
  9. Van der Tempel, L. Temperature solution for transient heat conduction in a thin bilayer between semi-infinite media in thermal effusivity measurement. J. Heat Transf. 2016, 138, 051301. [Google Scholar] [CrossRef]
  10. Zhou, L.; Bai, M.; Cui, W.; Lü, J. Theoretical solution of transient heat conduction problem in one-dimensional multilayer composite media. J. Thermophys. Heat Transf. 2017, 31, 483–488. [Google Scholar] [CrossRef]
  11. Cole, K.D.; Beck, J.V.; Woodbury, A.; de Monte, F. Intrinsic verification and a heat conduction database. Int. J. Therm. Sci. 2014, 78, 36–47. [Google Scholar] [CrossRef]
  12. De Monte, F.; Beck, J.V. X40B1T00, Exact Analytical Conduction Toolbox. October 2013, pp. 1–18. Available online: www.exact.unl.edu (accessed on 1 October 2020).
  13. Giacchetti, L.; de Monte, F. X60B1T00, Exact Analytical Conduction Toolbox. September 2015, pp. 1–72. Available online: www.exact.unl.edu (accessed on 1 October 2020).
  14. Haji-Sheikh, A.; Beck, J.V. An efficient method of computing eigenvalues in heat conduction. Numer. Heat Transf. Part B Fundam. 2000, 38, 133–156. [Google Scholar] [CrossRef]
  15. D’Alessandro, G.; de Monte, F. X42B10T00, Exact Analytical Conduction Toolbox. July 2016, pp. 1–73. Available online: www.exact.unl.edu (accessed on 1 October 2020).
  16. De Monte, F.; Beck, J.V.; Amos, D.E. A heat-flux based building block approach for solving heat conduction problems. Int. J. Heat Mass Transf. 2011, 54, 2789–2800. [Google Scholar] [CrossRef]
  17. De Monte, F.; Beck, J.V.; Amos, D.E. Solving two-dimensional Cartesian unsteady heat conduction problems for small values of the time. Int. J. Therm. Sci. 2012, 60, 106–113. [Google Scholar] [CrossRef]
  18. D’Alessandro, G.; de Monte, F. Sensitivity coefficients for thermal properties measurements using a boundary condition of the 4th kind. J. Phys. Conf. Ser. 2017, 796, 012011. [Google Scholar] [CrossRef]
  19. Arpaci, V.S. Conduction Heat Transfer; Addison-Wesley Publishing Company: Ann Arbor, MI, USA, 1966. [Google Scholar]
  20. D’Alessandro, G.; de Monte, F.; Amos, D. X62B10T00, Exact Analytical Conduction Toolbox. February 2017, pp. 1–95. Available online: www.exact.unl.edu (accessed on 1 October 2020).
  21. D’Alessandro, G.; de Monte, F. On the optimum experiment and heating times when estimating thermal properties through the plane source method. Heat Transf. Eng. 2021, 43. in press. [Google Scholar]
  22. De Monte, F.; Beck, J.V. X20B1T0, Exact Analytical Conduction Toolbox. January 2013, pp. 1–12. Available online: www.exact.unl.edu (accessed on 30 November 2020).
  23. De Monte, F.; Beck, J.V. X22B10T0, Exact Analytical Conduction Toolbox. February 2013, pp. 1–14. Available online: www.exact.unl.edu (accessed on 30 November 2020).
Figure 1. Schematic of the X40B1T11 problem.
Figure 1. Schematic of the X40B1T11 problem.
Energies 13 06484 g001
Figure 2. Three-layer structure: thin layer having 2Lf as thickness, in perfect contact with two semi-infinite bodies.
Figure 2. Three-layer structure: thin layer having 2Lf as thickness, in perfect contact with two semi-infinite bodies.
Energies 13 06484 g002
Figure 3. Schematic of the X42B10T11 problem.
Figure 3. Schematic of the X42B10T11 problem.
Energies 13 06484 g003
Figure 4. Three-layer structure: thin layer having 2Lf as thickness, in perfect contact with two slabs.
Figure 4. Three-layer structure: thin layer having 2Lf as thickness, in perfect contact with two slabs.
Energies 13 06484 g004
Figure 5. Dimensionless temperature plots for the X40B1T00 problem: (a) Temperature vs. x ˜ and t ˜ as a parameter; (b) Temperature vs. t ˜ and x ˜ as a parameter.
Figure 5. Dimensionless temperature plots for the X40B1T00 problem: (a) Temperature vs. x ˜ and t ˜ as a parameter; (b) Temperature vs. t ˜ and x ˜ as a parameter.
Energies 13 06484 g005
Figure 6. Dimensionless temperature distribution for the X40B1T00 problem for different values of the heat capacity ratio.
Figure 6. Dimensionless temperature distribution for the X40B1T00 problem for different values of the heat capacity ratio.
Energies 13 06484 g006
Figure 7. Dimensionless temperature plots for the X60B1T00 problem: (a) Temperature vs. x ˜ and t ˜ as a parameter; (b) Temperature vs. t ˜ and x ˜ as a parameter.
Figure 7. Dimensionless temperature plots for the X60B1T00 problem: (a) Temperature vs. x ˜ and t ˜ as a parameter; (b) Temperature vs. t ˜ and x ˜ as a parameter.
Energies 13 06484 g007
Figure 8. Dimensionless temperature distribution for the X60B1T00 problem for different values of the dimensionless contact resistance.
Figure 8. Dimensionless temperature distribution for the X60B1T00 problem for different values of the dimensionless contact resistance.
Energies 13 06484 g008
Figure 9. Dimensionless temperature plots for the X42B10T00 problem: (a) Temperature vs. x ˜ and t ˜ as a parameter; (b) Temperature vs. t ˜ and x ˜ as a parameter.
Figure 9. Dimensionless temperature plots for the X42B10T00 problem: (a) Temperature vs. x ˜ and t ˜ as a parameter; (b) Temperature vs. t ˜ and x ˜ as a parameter.
Energies 13 06484 g009
Figure 10. Dimensionless temperature distribution for the X42B10T00 problem for different values of the heat capacity ratio ( t ˜ = 0.2 ).
Figure 10. Dimensionless temperature distribution for the X42B10T00 problem for different values of the heat capacity ratio ( t ˜ = 0.2 ).
Energies 13 06484 g010
Figure 11. Dimensionless temperature plots for the X62B10T00 problem: (a) Temperature vs. x ˜ and t ˜ as a parameter; (b) Temperature vs. t ˜ and x ˜ as a parameter (P = 1 and R ˜ c = 0.25 ).
Figure 11. Dimensionless temperature plots for the X62B10T00 problem: (a) Temperature vs. x ˜ and t ˜ as a parameter; (b) Temperature vs. t ˜ and x ˜ as a parameter (P = 1 and R ˜ c = 0.25 ).
Energies 13 06484 g011
Figure 12. Dimensionless temperature distribution for the X62B10T00 problem for different values of the dimensionless contact resistance and P = 0.1.
Figure 12. Dimensionless temperature distribution for the X62B10T00 problem for different values of the dimensionless contact resistance and P = 0.1.
Energies 13 06484 g012
Table 1. Deviation time and quasi-steady times for different numerical accuracies A.
Table 1. Deviation time and quasi-steady times for different numerical accuracies A.
A t ˜ d   ( x ˜ , A ) t ˜ s   ( P , A )
x ˜ = 0 x ˜ = 0.25 x ˜ = 0.5 x ˜ = 1 P = 0.01 P = 0.1 P = 1 P = 10
20.20000.15310.11250.05002.23462.30382.78844.1702
50.08000.06130.04500.02005.03425.10345.58806.9698
100.04000.03060.02250.01009.70029.769410.254011.6358
150.02670.02040.01500.006714.366214.435414.920016.3018
Table 2. Quasi-steady times for several P and R ˜ c values and two different accuracies A.
Table 2. Quasi-steady times for several P and R ˜ c values and two different accuracies A.
P A = 2 A = 10
R ˜ c = 0.01 R ˜ c = 0.1 R ˜ c = 1 R ˜ c = 0.01 R ˜ c = 0.1 R ˜ c = 1
0.010.95200.95210.95403.77193.77263.7798
0.11.12541.14051.43964.45914.51905.7040
12.26812.56096.31388.986810.147025.0168
103.52494.153011.328113.966416.455344.8845
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

D’Alessandro, G.; de Monte, F. Multi-Layer Transient Heat Conduction Involving Perfectly-Conducting Solids. Energies 2020, 13, 6484. https://doi.org/10.3390/en13246484

AMA Style

D’Alessandro G, de Monte F. Multi-Layer Transient Heat Conduction Involving Perfectly-Conducting Solids. Energies. 2020; 13(24):6484. https://doi.org/10.3390/en13246484

Chicago/Turabian Style

D’Alessandro, Giampaolo, and Filippo de Monte. 2020. "Multi-Layer Transient Heat Conduction Involving Perfectly-Conducting Solids" Energies 13, no. 24: 6484. https://doi.org/10.3390/en13246484

APA Style

D’Alessandro, G., & de Monte, F. (2020). Multi-Layer Transient Heat Conduction Involving Perfectly-Conducting Solids. Energies, 13(24), 6484. https://doi.org/10.3390/en13246484

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop