Construction of Short-Time Heat Conduction Solutions in One-Dimensional Finite Rectangular Bodies

The symbolic expression of an “exact” analytical solution to a transient, 1D heat conduction problem of Fourier type for a finite body, denoted by $Te(x,t)$⁠, with $x∈[0,L]$ and $t≥0$⁠, is available only in an infinite series form, that is, consisting of an infinite number of terms. Though the solution to the above initial-boundary value problem is unique [1 (Chapter 1)], the symbolic expression of the related exact analytical solution is not unique. In fact, there exist two different symbolic forms of it in the specialized heat conduction literature [2–4]: (1) the short-time form coming from Laplace transform combined to an asymptotic series expansion and (2) the large-time form deriving from the application of the separation-of-variables method. Both short- and large-time expressions of the solution exhibit single summations consisting of an infinite number of terms.

As no calculator can account for an infinite number of terms, the exact analytical solution needs to be replaced with an approximate one consisting of a finite number of terms. Though this is well-known, only a few existing references address the required number of terms in the series evaluation for small and large values of the time. It seems that the mathematical expression is all the reader needs, when in fact, for engineering applications, considerable additional work is usually required to obtain reliable and accurate numerical values.

In an earlier paper by Beck et al. [5], the possibility of using small-time approximate solutions at the heated surface of some basic one-dimensional geometries was investigated and discussed. The treatment was, however, limited to the computation of only surface temperature for only constant-heat-flux boundary conditions. Also, for a plate, the approximate solution was consisting of two terms and valid for dimensionless times $t˜=αt/L2$ less than 0.2. But its accuracy in this time range was not given. A related paper restricted to 1D simple geometries was presented by Lavine and Bergman [6]. The small-time approximations derived by the authors using the Euler-Maclaurin summation formula for only surface temperature coincide with the 1D semi-infinite solutions for plane walls and, hence, it consists of only one term. It is suggested to use it for dimensionless times less than 0.2 but no accuracy is given, and only constant surface heat flux condition was investigated.

Another related paper was presented later by de-Monte et al. [7] for only rectangular bodies. But the proposed small-time approximate solutions are given for both temperature and heat flux and are valid for any location of the body and for any kind of boundary condition and its time variation. In addition, their time range of applicability is given as a function of location of evaluation and accuracy desired there. In detail, for a 1D finite rectangular body, temperature and heat flux can be assumed to be zero at a location $x$ with errors less than 10^-A, being A = 2, 3, …, 15, if there the time is less than the so-called penetration time $tp=x2/(10Aα)$⁠. It is defined as the time that it takes for a point of a body to be just affected by the “active” heated/cooled boundary in terms of temperature and heat flux. By “just affected” it means at the level of one part in 10^A with respect to the maximum temperature rise/fall that occurs at $x=0$ at the same time, i.e., $ΔTmax=|T(0,tp)−Tin|$⁠, where $Tin$ is the uniform initial temperature of the plate. In addition, temperature and heat flux can be computed at a location $x$ of the plate using the solutions of the related semi-infinite problem with errors less than 10^-A, if there the time is less than the so-called (first) deviation time $td(1)=(2L−x)2/(10Aα)≥tp$⁠. It is defined as the time that it takes for a point of a plane wall heated/cooled at $x=0$ to be affected (in terms of temperature and heat flux) by the presence of the homogeneous (inactive) boundary at $x=L$ at the level of one part in 10^A with respect to $ΔTmax$⁠.

A further paper on the same topic was proposed by McMasters et al. [8]. It extends the concept of penetration time given in Ref. [7] to cylindrical and spherical geometries. In addition, it is proven that the penetration time corresponding to a constant-temperature boundary condition, defined through $tp=x2/(10Aα)$ for rectangular 1D bodies [7], is the most conservative among the $T0(t/t0)n/2$ time-varying boundary conditions, with n = 0, 1,2, …, and $t0$ arbitrary reference time. By extension, the first kind boundary condition can also be considered the most conservative case for spherical and cylindrical geometries. Later Ostrogorsky and Mikic [9] extend the semi-infinite solutions to radial systems such as cylinders and spheres for short values of the time. Boundary conditions of first, second and third kind are considered. The time ranges of validity of the small-time approximations consisting of only one term and applicable at any location of the body are given. Also, their accuracy is provided in a graphical form through a relative error.

The objective of this paper is to construct very simple analytical solutions for 1D finite rectangular bodies applicable at any location and at the level of accuracy desired (i.e., one part in $10A$ with respect to $ΔTmax$ defined before) for small values of the time but larger than the (first) deviation time $td(1)$ defined in the previous subsection. The extension of the time range of applicability of the solution requires increasing its number of terms consisting currently of only one single semi-infinite term for times less than or equal to $td(1)$ [7]. For this purpose, a new short-time solution will be constructed superimposing two appropriate semi-infinite terms and, hence, will be valid for longer times till $t=td(2)>td(1)$⁠. An algebraic conservative expression is also derived for $td(2)$ (termed as “second deviation time”) depending on both the location of evaluation and accuracy desired.

Therefore, the time interval extension stated before is the prime objective of the current work as it allows matching better from a computational viewpoint the short- and large-time forms of the solution, so facilitating time partitioning. For this reason, the second deviation time may also be considered as a partitioning time. A first paper dealing with the analysis and implementation of the time partitioning method was proposed by Yen et al. [10] where, however, an algebraic expression was not given. In fact, the optimum partitioning time was estimated only a posteriori through several numerical experiments. Then, it was used for intrinsic verification purposes of numerical methods in transient heat conduction [11] and for improving convergence of infinite series involved in transient temperature solutions [4 (Sec. 5.3)].

The construction of the two-term short-time solution may be obtained by using: (1) some adequate physical considerations (due to thermal symmetries with respect to the middle plane of a slab) to simulate homogeneous boundary conditions of the first and second kinds on the back side of the slab (inactive surface); (2) the principle of superposition for linear problems [12 (Chapter 3)]; and (3) the concepts of penetration time and (first) deviation time that allow the well-known solutions of a semi-infinite body to be utilized as basic “building blocks” for constructing the short-time solution in terms of temperature and heat flux. Before proceeding to apply this methodology, the concept of “computational” analytical solution versus “exact” analytical solution and related relative error will be defined in Sec. 2.

The method here proposed applies to transient problems with active boundary conditions of the first, second, and third kinds on one side of the slab, and inactive boundary conditions of only the first and second kinds on the other side. In particular, the proposed method is used to find short-time solutions applicable at different times and different locations within the body with the accuracy desired for six different flat plates cases (Secs. 3–8). They are denoted by X11, X12, X21, X22, X31, and X32 in the heat conduction numbering system devised by Cole et al. [4 (Chapter 2)], where 1, 2, and 3 indicate the boundary conditions of first, second, and third kinds, respectively. See also Ref. [12 (Chapter 2 and Appendix A)] by Woodbury et al. Then, two examples on the maximum number of terms required to obtain a temperature solution with three different accuracies matching the two-term short-time form here derived and the well-known large-time form are given in Sec. 9. Some plots and tables are also provided.

The method is illustrated for constant boundary conditions in Secs. 3–8 but can be extended to any time-dependent boundary function applied at the active boundary surface, as shown in Sec. 10.

where $E=E(x,t,M)$ is hence the relative error.

where $T0(x,t)$ and $Tm(x,t)$ depend on (1) the form chosen for the analytical solution (short-time form or large-time form), (2) boundary condition type, and (3) boundary function prescribed.

Equation (2) states that, when $M→∞$⁠, $Tc(x,t,M)→Te(x,t)$⁠. Also, it shows that the numerator of Eq. (1) is the absolute value of the “tail” of the summation. This tail consisting of an infinite number of terms may be calculated conservatively by using the Euler-Maclaurin formula [13 (Eq. (4.14.1))]. An example of its use is given in Ref. [4 (Sec. 5.1.3)] for large-time solutions. The above formula can also be applied to the finite summation appearing at the denominator of Eq. (2).

By setting $E=10−A$ in Eq. (1), where A = 2, 3, …, 15 is the counting integer for the accuracy desired, and solving the resulting equation, the maximum number of terms $MA$ to get an accuracy of one part in $10A$ (with respect to the maximum temperature change $ΔTmax=|Tc(0,t,M)−Tin|$⁠) may be calculated as a function of time $t$⁠, location $x$⁠, and accuracy A. In detail, A = 2 (accuracy of 1%) is for visual comparison and is acceptable in many engineering applications, while A = 10 or even 15 is mainly for verification purposes of large numerical codes [14–16]. Note that the floating-point relative accuracy in double precision within MATLAB software is of 2.2204 × 10⁻¹⁶. When using A = 15, the “computational” analytical solution may be considered in every practical respect as the “exact” one.

where $Tc(0,t,MA)=Te(0,t)$ when the flat plate is subject to a boundary condition of the first kind.

Similarly, an error between computational and exact analytical heat flux solutions can be defined.

Consider a slab of thickness L with temperature-independent properties, initially at zero temperature (without any loss of generality), subject to a step change in temperature (⁠$0→T0$⁠) at x = 0 and thermally insulated at x = L, as shown in Fig. 1, where P denotes a generic location of abscissa x. By using the numbering system devised by Cole et al. [4 (Chapter 2)], this problem may be denoted by X12B10T0, where “1” in B10 denotes a constant boundary condition at x = 0, while “0” indicates a homogeneous boundary condition at x = L. Also, “T0” refers to a zero initial temperature within the plate.

Very simple solutions consisting of only two terms can be constructed and, hence, used to compute temperature and heat flux at any location of the slab at early times with the accuracy desired in place of the exact analytical solutions consisting of infinite terms, as shown in next subsections.

with a relative error (as defined by Eq. (1) for temperature) less than 10^−A.

which are the well-known solutions of the X10B1T0 problem, being $erfc(.)$ the complementary error function [4 (Appendix E)].

The slab of Fig. 1 is equivalent to a slab 2 L thick and boundary condition like that at x = 0 imposed at x = 2 L, i.e., X12B10T0 = X11B11T0, as shown in Fig. 3. In fact, if the slab is subject to a second kind of homogeneous boundary condition at x = L (Fig. 1), it can be simulated by a second layer in perfect thermal contact with the first one and having its same thickness, and subject to a nonzero boundary term at x = 2 L that is equal to that at x = 0, as shown in Fig. 3.

where X11B10T0 and X11B01T0 are depicted in Figs. 4(a) and 4(b), respectively.

where 2 L has been replaced by 4 L, the slab may be approximated to a semi-infinite rectangular body (⁠$x≥0$⁠), i.e., X11B10T0 ≈ X10B1T0 (see Fig. 5(a)), with relative errors E less than 10^−A.

which are the solutions of the X10B1T0 problem of Fig. 5(a).

where 2 L − x has been replaced by 2 L + x, the slab may be approximate to a semi-infinite rectangular body (⁠$x≤2L$⁠), i.e., X11B01T0 ≈ X01B1T0 (see Fig. 5(b)), with relative errors E less than 10^−A.

which are the solutions of the X01B1T0 problem of Fig. 5(b) related to the ones of the companion X10B1T0 case of Fig. 5(a).

Also, note that $td(1a)$ defined by Eq. (9) is always greater than $td(1b)$ given by Eq. (11) with the only exception of the inactive back side x = L of the slab of Fig. 1, where the above two deviation times are the same. As $td(1b)=min{td(1a),td(1b)}$⁠, for conservative purposes it may be taken as second deviation time, that is, $td(2)=td(1b)$⁠.

Consider a slab of thickness L with temperature-independent properties, initially at zero temperature, subject to a step change in temperature (⁠$0→T0$⁠) at x = 0 and kept at zero temperature at x = L, as shown in Fig. 6 where P denotes a generic location of abscissa x. By using the numbering system [4 (Chapter 2)], this problem may be denoted by X11B10T0.

Very simple solutions consisting of only two terms can be constructed and, hence, used to compute temperature and heat flux at any location of the slab at early times with the accuracy desired in place of the exact analytical solutions consisting of infinite terms, as shown in next subsections.

The short-time computational solutions to the current problem, with a relative error $E≤10−A$⁠, being A = 2, 3, …, 15, are [7]

• $0≤t≤tp$⁠, where $tp=x2/(10Aα)$

  $TX11B10T0(x,t)|0≤x≤L≈0$

  (14a)

  $qX11B10T0(x,t)|0≤x≤L≈0$

  (14b)
• $tp<t≤td(1)=(2L−x)2/(10Aα)$⁠:

  $TX11B10T0(x,t)|0≤x≤L≈T0erfc(x2αt)$

  (15a)

  $qX11B10T0(x,t)|0≤x≤L≈kT0παte−x24αt$

  (15b)

The slab of Fig. 6 is equivalent to a slab 2 L thick and boundary condition opposite to that at x = 0 imposed at x = 2 L, i.e., X11B10T0 = X11B11T0, as shown in Fig. 7. In fact, if the slab is subject to a first kind of homogeneous boundary condition at x = L (Fig. 6), it can be simulated by a second layer in perfect thermal contact with the first one and having its same thickness, and subject to a nonzero boundary term at x = 2 L that is opposite to that at x = 0, as shown in Fig. 7.

where X11B10T0 and X11B01T0 are depicted in Figs. 4(a) and 4(b), respectively.

Consider now the X11B10T0 problem of Fig. 4(a). If $0≤t≤td(1a)$⁠, where $td(1a)$ is the deviation time given by Eq. (9), the slab may be approximated to a semi-infinite rectangular body (⁠$x≥0$⁠), i.e., X11B10T0 ≈ X10B1T0 (see Fig. 5(a)), with relative errors E less than 10^−A.

Therefore, for $0≤x≤2L$ and $0≤t≤td(1a)$⁠, temperature and heat flux may be calculated at any point P of abscissa x using Eqs. (10a) and (10b), which are the solutions of the X10B1T0 problem of Fig. 5(a).

Consider now the X11B01T0 problem of Fig. 4(b). If $0≤t≤td(1b)$⁠, where $td(1b)$ is the deviation time given by Eq. (11), the slab may be approximate to a semi-infinite rectangular body (⁠$x≤2L$⁠), i.e., X11B01T0 ≈ X01B1T0 (see Fig. 5(b)), with relative errors E less than 10^−A.

Therefore, by changing $x→(2L−x)$ in Eqs. (10a) and (10b), for $0≤x≤2L$ and $0≤t≤td(1b)$ the temperature and heat flux may be calculated through Eqs. (12a) and (12b), respectively, which are the solutions of the X01B1T0 problem of Fig. 5(b).

Consider a slab of thickness L with temperature-independent properties, initially at zero temperature, subject to a step change in heat flux (⁠$0→q0$⁠) at x = 0 and thermally insulated at x = L, as shown in Fig. 1, where $T0$ has to be replaced with $q0$⁠. By using the numbering system, this problem may be denoted by X22B10T0.

The short-time computational solutions to the current problem for $0≤x≤L$ and for a relative error $E≤10−A$⁠, with A = 2, 3, …, 15, are

• $0≤t≤tp$⁠, where $tp=x2/(10Aα)$

  $TX22B10T0(x,t)|0≤x≤L≈0$

  (20a)

  $qX22B10T0(x,t)|0≤x≤L≈0$

  (20b)
• $tp<t≤td(1)$⁠, where $td(1)=(2L−x)2/(10Aα)$

  $TX22B10T0(x,t)|0≤x≤L≈TX20B1T0(x,t)|x≥0$

  (21a)

  $qX22B10T0(x,t)|0≤x≤L≈qX20B1T0(x,t)|x≥0$

  (21b)
• $td(1)<t≤td(2)$⁠, where $td(2)=(2L+x)2/(10Aα)$

  $TX22B10T0(x,t)|0≤x≤L≈TX20B1T0(x,t)|x≥0+TX20B1T0(2L−x,t)|x≥0$

  (22a)

  $qX22B10T0(x,t)|0≤x≤L≈qX20B1T0(x,t)|x≥0−qX20B1T0(2L−x,t)|x≥0$

  (22b)

where the $ierfc(.)$function is the complementary error function integral [4 (Appendix E)].

Consider a slab of thickness L with temperature-independent properties, initially at zero temperature, subject to a step change in heat flux (⁠$0→q0$⁠) at x = 0 and kept at zero temperature at x = L, as shown in Fig. 6, where $T0$ has to be replaced with $q0$⁠. By using the numbering system, this problem may be denoted by X21B10T0.

The short-time computational solutions to the current problem for $0≤x≤L$ and for a relative error $E≤10−A$⁠, with A = 2, 3, …, 15, are (see Sec. 4.1)

• $0≤t≤tp$⁠, where $tp=x2/(10Aα)$

  $TX21B10T0(x,t)|0≤x≤L≈0$

  (24a)

  $qX21B10T0(x,t)|0≤x≤L≈0$

  (24b)
• $tp<t≤td(1)$⁠, where $td(1)=(2L−x)2/(10Aα)$

  $TX21B10T0(x,t)|0≤x≤L≈TX20B1T0(x,t)|x≥0$

  (25a)

  $qX21B10T0(x,t)|0≤x≤L≈qX20B1T0(x,t)|x≥0$

  (25b)
• $td(1)<t≤td(2)$⁠, where $td(2)=(2L+x)2/(10Aα)$

  $TX21B10T0(x,t)|0≤x≤L≈TX20B1T0(x,t)|x≥0−TX20B1T0(2L−x,t)|x≥0$

  (26a)

  $qX21B10T0(x,t)|0≤x≤L≈qX20B1T0(x,t)|x≥0+qX20B1T0(2L−x,t)|x≥0$

  (26b)

where $TX20B1T0(x,t)|x≥0$ and $qX20B1T0(x,t)|x≥0$ are given by Eqs. (23a) and (23b), respectively.

Consider a slab of thickness L with temperature-independent properties, initially at zero temperature, subject to a step change in fluid temperature (⁠$0→T∞,0$⁠) at x = 0 and thermally insulated at x = L, as shown in Fig. 1, where $T0$ has to be replaced with $T∞,0$⁠. Also, the heat exchange by convection at that boundary is characterized by a coefficient $h0$⁠. By using the numbering system, this problem may be denoted by X32B10T0.

The short-time computational solutions to the current problem for $0≤x≤L$ and for a relative error $E≤10−A$⁠, with A = 2, 3, …, 15, are

• $0≤t≤tp$⁠, where $tp=x2/(10Aα)$⁠:

  $TX32B10T0(x,t)|0≤x≤L≈0$

  (27a)

  $qX32B10T0(x,t)|0≤x≤L≈0$

  (27b)
• $tp<t≤td(1)$⁠, where $td(1)=(2L−x)2/(10Aα)$⁠:

  $TX32B10T0(x,t)|0≤x≤L≈TX30B1T0(x,t)|x≥0$

  (28a)

  $qX32B10T0(x,t)|0≤x≤L≈qX30B1T0(x,t)|x≥0$

  (28b)
• $td(1)<t≤td(2)$⁠, where $td(2)=(2L+x)2/(10Aα)$

  $TX32B10T0(x,t)|0≤x≤L≈TX30B1T0(x,t)|x≥0+TX30B1T0(2L−x,t)|x≥0$

  (29a)

  $qX32B10T0(x,t)|0≤x≤L≈qX30B1T0(x,t)|x≥0−qX30B1T0(2L−x,t)|x≥0$

  (29b)

The $erfcx(z)=ez2erfc(z)$ is the scaled complementary error function available in Matlab ambient. It is recommended for large values of x, t, and h₀ as difficulties might arise during the numerical evaluation of Eq. (30c) in its former expression because of the positive argument of the exponential term (depending on the variables listed before). On the contrary, the scaled complementary error function (appearing in the latter expression) approaches zero for large values of x, t and h₀. A proof of this statement may be derived using the asymptotic expansion of [4 (Eq. (E.6))].

Consider a slab of thickness L with temperature-independent properties, initially at zero temperature, subject to a step change in fluid temperature (⁠$0→T∞,0$⁠) at x = 0 and kept at zero temperature at x = L, as shown in Fig. 6, where $T0$ has to be replaced with $T∞,0$⁠. Also, the heat exchange by convection is characterized by $h0$⁠. By using the numbering system, this problem may be denoted by X31B10T0.

The short-time computational solutions to the current problem for $0≤x≤L$ and for a relative error $E≤10−A$⁠, with A = 2, 3, …, 15, are

• $0≤t≤tp$⁠, where $tp=x2/(10Aα)$

  $TX31B10T0(x,t)|0≤x≤L≈0$

  (31a)

  $qX31B10T0(x,t)|0≤x≤L≈0$

  (31b)
• $tp<t≤td(1)$⁠, where $td(1)=(2L−x)2/(10Aα)$

  $TX31B10T0(x,t)|0≤x≤L≈TX30B1T0(x,t)|x≥0$

  (32a)

  $qX31B10T0(x,t)|0≤x≤L≈qX30B1T0(x,t)|x≥0$

  (32b)
• $td(1)<t≤td(2)$⁠, where $td(2)=(2L+x)2/(10Aα)$⁠:

  $TX31B10T0(x,t)|0≤x≤L≈TX30B1T0(x,t)|x≥0−TX30B1T0(2L−x,t)|x≥0$

  (33a)

  $qX31B10T0(x,t)|0≤x≤L≈qX30B1T0(x,t)|x≥0+qX30B1T0(2L−x,t)|x≥0$

  (33b)

where $TX30B1T0(x,t)|x≥0$ and $qX30B1T0(x,t)|x≥0$ are given by Eqs. (30a) and (30b), respectively, along with Eq. (30c).

The penetration time and the first and second deviation times defined in the previous sections are given in dimensionless form by Table 1 for different values of the accuracy A and at two different locations, namely, the front side $x˜=0$ of the slab and its back side $x˜=1$⁠, with $x˜=x/L$⁠. In detail, $t˜p=x˜2/(10A)$⁠, $t˜d(1)=(2−x˜)2/(10A)$⁠, and $t˜d(2)=(2+x˜)2/(10A)$⁠, where the dimensionless time may be defined as $t˜=αt/L2$⁠.