- Email: [email protected]
Dynamic response of a non-homogeneous 1D slab under periodic thermal excitation
International Journal of Heat and Mass Transfer 50 (2007) 3943–3948 www.elsevier.com/locate/ijhmt
Dynamic response of a non-homogeneous 1D slab under periodic thermal excitation G.E. Cossali * Facolta` di Ingegneria, Universita` di Bergamo, via Marconi 6, 24044 Dalmine (BG), Italy Received 6 June 2006 Available online 30 March 2007
Abstract The periodic heat conduction in a non-homogeneous slab is studied and an approximate analytic solution, valid for weakly non homogeneous slabs, is found. An exact solution for a well defined non-homogeneity is also found and compared to the approximate solution to estimate the accuracy of the result. The solution allows to extend the method of thermal quadrupoles to weakly non-homogeneous 1D-slabs. Ó 2007 Elsevier Ltd. All rights reserved.
1. Introduction Transient heat conduction in homogeneous materials is a deeply analysed topic, also for the vast consequences in many applied fields, like energy storage [1], laser heating, thermal behaviour of building materials [2], but also for measurement techniques like microcalorimetry [3] or transient techniques for measuring convective heat transfer coefficient [4,5]. Many different techniques were developed to solve transient one-dimensional problems also for the more complex case of composite walls [6]. The case of transient heating of a non-homogeneous materials is more complex and numerical approaches were used to find solutions for a certain class of materials [7,8]. Periodic heat conduction is a relatively less investigated topic, although the implications for many applied fields are significant, for example for building walls and materials [9,10] or again for techniques to measure conductivity and heat capacity like a.c. calorimetry [11] or 3-x methods [12] and for local convective heat transfer coefficient measurements [13]. The periodic conduction in a 1D homogeneous slab is a relatively simple problem but recently many studies, dealing
with it, can be found in the literature [16–19], particularly related to the hyperbolic version of the heat equation, as the wave like behaviour can be better evidenced in a simple geometry. However, the extension of the transient conduction analysis to non-homogeneous material is not straightforward, and fewer works can be found on this subject [14,15], whereas periodic conduction in non-homogeneous material does not seems to have been sufficiently considered, although non-homogeneous material are becoming everyday more widespread in many engineering applications. The paper presents a solution of the periodic conduction in a non-homogeneous 1D slab under the condition of weak non-homogeneity (i.e. limited variation of the material properties with position) and an assessment of the accuracy of the method by comparison with an exact solution. The solution allows a relatively simple extension of the thermal quadrupole method to this class of materials. 2. Basic equations Consider a 1D slab made of a non-homogeneous material, the energy equation can then be written as qc
*
Tel.: +39 03520 52309; fax: +39 03556 2779. E-mail address: [email protected]
0017-9310/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2007.01.038
oT oq ¼ ot ox
and making use of Fourier law: q ¼ k oT we obtain: ox
ð1Þ
3944
d2
G.E. Cossali / International Journal of Heat and Mass Transfer 50 (2007) 3943–3948
oT oT o2 T ¼ cðnÞ þ aðnÞ 2 ot on on
ð2Þ
ok ; and aðnÞ ¼ qck . The temperature where n ¼ dx ; cðnÞ ¼ qc1 on and heat flux fields can be split into time average (T a ; qa ) and fluctuating ( Te , ~ q) components. Introducing the Fourier transform of the latter: Z þ1 T ðx; tÞ ¼ T a ðxÞ þ Sðx; xÞeixt dx; 1 Z þ1 Qðx; xÞeixt dx qðx; tÞ ¼ qa ðxÞ þ 1
Eq. (2) becomes: ixd2 S ¼ cðnÞS 0 þ aðnÞS 00
ð3Þ
0 ¼ cðnÞT 0a þ aðnÞT 00a
ð4Þ
where apex means derivation respect to n. The case of homogeneous slab is found setting c ¼ 0 and a ¼ a0 ¼ const: and Eq. (3) becomes: S 00H b2 S H ¼ 0 qffiffiffiffiffiffiffi 2 with b ¼ ixd , the solution is easily found under the more a0 convenient form: S H ¼ S H;0 cosh ðbnÞ þ
S 0H;0 sinh ðbnÞ b
ð5Þ
with S H;0 ¼ S H ð0Þ and S 0H;0 ¼ S 0H ð0Þ. The general solution of Eq. (4) is easily found by double integration: Z n R f cðvÞ dv Ta ¼ A þ B e 0 aðvÞ df: 0
mðnÞ ¼ S H0 m1 ðnÞ þ
S 0H0 m2 ðnÞ b
ð7Þ
cðnÞ aðnÞ 2 m1 ðnÞ ¼ b sinh ðbnÞ þ b cosh ðbnÞ ; a0 a0 cðnÞ aðnÞ 2 m2 ðnÞ ¼ b coshðbnÞ þ b sinhðbnÞ a0 a0 are known functions. The solution of (6) is then: Z 1 1 n s ¼ s0 coshðbnÞ þ s00 sinhðbnÞ mðzÞ sinh½bðz nÞdz b b 0 and choosing sð0Þ ¼ 0, s0 ð0Þ ¼ 0, we obtain: Z 1 n S ¼ S 0 cosh ðbnÞ m1 ðzÞ sinh½bðz nÞdz b 0 Z S0 1 n þ 0 sinh ðbnÞ m2 ðzÞ sinh½bðz nÞdz b 0 b
ð8Þ
with S 0 ¼ Sð0Þ and S 00 ¼ S 0 ð0Þ. The linearisation that led to Eq. (6) introduces approximation whose entity depends on the magnitude of the slab thermal properties variations. To assess the validity of the approximation a comparison to an exact solution of Eq. (3), for a particular distribution of the slab properties, is needed. 4. An exact solution Consider the following particular form for the dependence on position of the material properties: qc ¼ const; k ¼ Kðn0 nÞ
2
ð9Þ
Eq. (3) becomes qcxd2 S¼0 K and after the change of qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi variable: g ¼ ðn n0 Þ, and defin4xqcd2 1 ing k0 ¼ 2 1 þ i K , an Euler equation is promptly recovered: 1 2 0 2 00 k0 S ¼ 0 g S þ 2gS þ ð10Þ 4 2
ðn0 nÞ S 00 þ 2ðn n0 ÞS 0 i
3. The weakly non-homogeneous slab Consider a material characterised by variations of its properties with position that are small compared to the average value, i.e.: k ¼ k 0 þ jðnÞ;
jjðnÞj jk 0 j
with
0
0
jj ðnÞj ¼ jk ðnÞj jkðnÞj a ¼ a0 þ aðnÞ;
with
jaðnÞj ja0 j 0
from which it follows: j ac j ¼ j kk j 1. Suppose now that the solution of Eq. (3) can be written as a perturbation of the solution SH valid for homogeneous materials (Eq. (5)), i.e. S ¼ SH þ s
with
s ¼ oðS H Þ
then substituting into Eq. (3) and preserving only the terms up to the first order: s00
ixd 2 cðnÞ 0 aðnÞ 00 s¼ S S ¼ mðnÞ a0 H a0 H a0
where
ð6Þ
where now the apex means derivation respect to g. The solution of this equation can be written in the convenient form: 1 1 n0 n 2 1 S¼ S 0 sinhðzÞ þ k0 coshðzÞ S 00 n0 sinhðzÞ k0 n0 2 ð11Þ with
n n z ¼ k0 ln 0 n0
and again: Sð0Þ ¼ S 0 , S 0 ð0Þ ¼ S 00 . The exact analytic solution can now be used to assess the validity of the approximate solution (8).
G.E. Cossali / International Journal of Heat and Mass Transfer 50 (2007) 3943–3948
(
5. Assessment of the perturbative solution F 1;ap Consider again the material described by (9) with the further condition n0 1, then, after defining: k 0 ¼ Kn20 we obtain:
F 2;ap
) ! N 1 ðnÞ M 1 ðnÞ ¼ 1þ 2 sinh ðbnÞ cosh ðbnÞ þ n0 n20 ( ) " # 1 M 2 ðnÞ N 2 ðnÞ sinh ðbnÞ 1 þ ¼ þ 2 coshðbnÞ b n20 n0
jðnÞ ¼ k k 0 ¼ Kn2 2Kn0 n k 0 ðxÞ ¼ 2Kðn n0 Þ
ð13Þ with the functions:
and the condition: jk 0 j jkj yields: 2 jðn n0 Þj that is obviously verified as 0 6 n 6 1 and n0 1. Moreover a0 ¼
k 0 Kn20 ¼ ; qc qc
cðnÞ 2ðn n0 Þ ¼ ; a0 n20
3945
aðnÞ n2 2n0 n ¼ a0 n20
The weakly non-homogeneous approximation holds and the functions m1 and m2 (see Eq. (7)) now become: ( ) 2ðn n0 Þ n2 2n0 n 2 m1 ðnÞ ¼ b sinh ðbnÞ þ b coshðbnÞ ; n20 n20 ( ) 2ðn n0 Þ n2 2n0 n 2 b cosh ðbnÞ þ b sinhðbnÞ m2 ðnÞ ¼ n20 n20 qffiffiffiffiffiffiffiffi 2 with b ¼ i xd and the following integrals: a0 Z n m1 ðzÞ sinh ðb½z xÞ dz I 1 ðnÞ ¼
1 1 n0 n 2 1 n0 n F 1;ex ¼ sinh k0 ln k0 n0 2 n0 n0 n þk0 cosh k0 ln ð14Þ n0 1 n0 n 2 n0 n0 n F 2;ex ¼ sinh k0 ln n0 n0 k0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi respectively, with now k0 ¼ 12 1 þ 4b2 n20 . Fig. 1 reports a comparison for a fixed value of x, having introduced the non-homogeneity magnitude e defined as e¼
k max k min jð1 2n0 Þj ¼ k0 n20
0
I 2 ðnÞ ¼
Z
n
m2 ðzÞ sinh ðb½z xÞ dz
0
must be evaluated to find the solution in an explicit form. In Appendix, the integrals are calculated, finding: b fM 1 sinhðbnÞ þ N 1 cosh ðbnÞg n20 b I 2 ðnÞ ¼ 2 fM 2 sinhðbnÞ þ N 2 cosh ðbnÞg n0 I 1 ðnÞ ¼
with n0 1 n b nn 1 þ n 0 bn2 n3 ; N 1 ¼ 0 n2 6 2 4 2b 4b 2 1 n0 1 2 1 n0 b 2 b 3 n n n M2 ¼ 2 n n ; N2 ¼ 4 4b 2 6 2 4b M1 ¼
The approximate solution then becomes: ( ) ! N 1 ðnÞ M 1 ðnÞ S ¼ S H0 1þ sinhðbnÞ cosh ðbnÞ þ n20 n20 ( ) " # S 0H0 M 2 ðnÞ N 2 ðnÞ sinh ðbnÞ 1 þ þ þ 2 cosh ðbnÞ b n20 n0 ð12Þ The exact (11) and approximate (12) solutions can now be compared by comparing the functions:
Fig. 1. Relative discrepancies between the approximated and the exact values of: (a) the function F1 and (b) the function F2 versus the position along the slab, for x ¼ 1.
3946
G.E. Cossali / International Journal of Heat and Mass Transfer 50 (2007) 3943–3948
It appears that the maximum discrepancy occurs at n ¼ 1, thus it is interesting to compare the two solutions (by comparing the functions (13) and (14)) at the fixed position n ¼ 1 for different values of the exciting frequency. To this end, observe that by introducing the parameter ~e ¼ n1 0 ¼ e þ oðeÞ, the coefficients of the hyperbolic functions in F j;a evaluated for n ¼ 1 can be expressed as polynomials of second order in ~e: ~e ~e2 N1 ¼ 1 2 4 n20 ~e 1 ~e2 1 b M1 b þ C 12;a ¼ 2 ¼ 2 2b 3 2 b n0 2 ~e ~e 1 N2 1 1 C 21;a ¼ ¼ þ b n20 2 2 2b2 3 ! ~e ~e2 1 1 M2 1 1þ 2 ¼ 1 C 22;a ¼ þ1 b b 2 4 b2 n0
Moreover, defining: u ¼ k0 ln n0n1 b the functions F j;e 0 can be written (for n ¼ 1) as 1 1 sinh ðuÞ þ cosh ðuÞ cosh ðbÞ F 1;ex ¼ ð1 ~eÞ 2 2k0 1 1 2 þð1 ~eÞ coshðuÞ þ sinh ðuÞ sinh ðbÞ 2k0 C 11;a ¼ 1 þ
12
F 2;ex
ð1 ~eÞ ¼ ~ek0
sinh ðuÞ cosh ðbÞ 1
ð1 ~eÞ 2 þ cosh ðuÞ sinh ðbÞ ~ek0
!
Fig. 2. Relative discrepancies between the approximated and the exact values of: (a) the function F1 and (b) the function F2 versus nondimensional frequency (b) for n ¼ 1.
and the corresponding coefficients of the hyperbolic functions are (up to Second order in ~e): ~e ~e2 C 11;e ¼ 1 þ ð1 þ b2 Þ þ oð~e2 Þ 8 2 ~e 1 ~e2 7b 1 b þ C 12;e ¼ þ oð~e2 Þ 2 b 4 3 2b ~e ~e2 7b2 1 C 21;e ¼ þ þ þ oð~e2 Þ 2 4 3b2 2b2 ! ~e 1 2 b4 þ 3b2 1 1 2 1 þ ~e þ oð~e Þ C 22;e ¼ b 2 8 b2
ð15Þ 6. The thermal quadrupole approach for the non-homogeneous slab
a comparison with those relative to the approximate solution shows that the error is of order ~e2 , i.e. of order e2 ~e2 ð3 þ b2 Þ þ oð~e2 Þ 8 3~e2 1 b þ oð~e2 Þ ¼ 4 2b ~e2 1 ¼ 3 þ oð~e2 Þ 4 2b2 ~e2 b2 þ 4 ¼ þ oð~e2 Þ 2b 4
DC 11 ¼ C 11;a C 11;e ¼ DC 12 ¼ C 12;a C 12;e DC 21 ¼ C 21;a C 21;e DC 22 ¼ C 22;a C 22;e
Fig. 2 shows the comparison of the two solutions, the discrepancy between the approximate and the exact one remains limited also for relatively large (0.2) values of e. However, relative errors (DC jk =C jk;e ) are expected to increase for very large and very small values of the frequency x, as it stems from Eqs. (15) and (16).
ð16Þ
The thermal quadrupole approach [20] gives a compact and elegant way to solve transient problems for 1D homogeneous slabs and the application to multilayered slabs is straightforward. The previous results allow to extend the thermal quadrupole approach also to 1D non-homogeneous materials. Defining S Wðx; xÞ ¼ Q the following relation, that is the basis of the thermal quadrupole method, holds for a homogeneous slab: Wðd; xÞ ¼ Mh Wð0; xÞ where
G.E. Cossali / International Journal of Heat and Mass Transfer 50 (2007) 3943–3948
" Mh ¼
cosh ðbnÞ
1 sinh ðbnÞ kb
kb sinh ðbnÞ
cosh ðbnÞ
#
and using partial integration, the following recurrence relation is found:
For the non-homogeneous slab a similar equation can now be written, in fact: Z 1 n S ¼ S 0 cosh ðbnÞ m1 ðzÞ sinh½bðz nÞ dz b 0 Z Q0 1 n sinh ðbnÞ m2 ðzÞ sinh½bðz nÞ dz b 0 kb Z n Q ¼ S 0 kb sinh ðbnÞ k m1 ðzÞ cosh½bðz nÞ dz þ Q0 cosh ðbnÞ þ
1 b
Z
0 n
m2 ðzÞ cosh½bðz nÞ dz
"
coshðbnÞ
I n ¼
e2bx n n I x 2b n1 2b
where: Z I ¼ 0
x
expð2bzÞ dz ¼
0
e2bx 1 2b
and the closed form: I n ¼
nk n e2bx X ð1Þ n! n! nþ1 xk þ ð1Þ nþ1 2b k¼0 k!ð2bÞnk ð2bÞ
ð17Þ
0
and then defining: Mnh ¼
3947
kb1 sinhðbnÞ
#
kb sinhðbnÞ coshðbnÞ 2 Rn b1 0 m1 ðzÞ sinh½bðz nÞdz þ4 Rn k 0 m1 ðzÞ cosh½bðz nÞdz " 1 # b I 1 kb12 I 2 ¼ Mh þ kI 3 b1 I 4
1 kb2 1 b
Rn
Rn 0
0
m2 ðzÞ sinh½bðz nÞdz
m2 ðzÞ cosh½bðz nÞdz
3 5
is promptly found. Consider now the four integrals: Z x n I ss ¼ zn sinh ðbzÞ sinh ðbðz xÞÞ dz Z0 x n I cs ¼ zn cosh ðbzÞ sinh ðbðz xÞÞ dz 0 Z x I sc ¼ zn sinh ðbzÞ cosh ðbðz xÞÞ dz 0 Z x I cc ¼ zn cosh ðbzÞ cosh ðbðz xÞÞ dz 0
the following equation is obtained: Wðd; xÞ ¼ Mnh Wð0; xÞ
a simple calculation gives:
If the functions cðnÞ and aðnÞ can be expressed as polynomial in n, then the integrals Ij can be written in a closed form (see Appendix). These results allow to extend the use of the thermal quadrupole method to 1-d weakly non-homogeneous slabs.
xnþ1 xnþ1 coshðbxÞ; I ncs ¼ A sinhðbxÞ; 2ðn þ 1Þ 2ðn þ 1Þ xnþ1 xnþ1 I nsc ¼ A þ sinhðbxÞ; I ncc ¼ Aþ þ coshðbxÞ 2ðn þ 1Þ 2ðn þ 1Þ
7. Conclusions The thermal response of a weakly non-homogeneous slab to periodic thermal excitation was analysed through a Fourier transform approach and an analytic solution of the transformed problem, valid in the limit of small nonhomogeneity, was found. The solution can be written in closed form whenever the dependence of the wall thermal characteristics on the position can be represented in polynomial form. To assess the validity of the approximation, a comparison with an exact solution of the transformed problem is reported, showing that accuracy is of second order in the non-homogeneity parameter e. The solution allows also to extend the method of thermal quadrupoles to the weakly non-homogeneous slab thus extending the capability of the method. Appendix Defining the integrals: Z x In ¼ zn expð2bzÞ dz 0
I nss ¼ Aþ
ð18Þ where 1 bx A ¼ ½ebx I þ n e In 4 and using Eq. (17), after few manipulations, yields: 0 1 ! n n1 2 2 X 1 X 1 Aþ ¼ xn2s an;n2s sinhðbxÞ þ @ xn2sþ1 an;n2sþ1 A coshðbxÞ 4b s¼0 4b s¼0 sinhðbxÞ n even n! þ 2ð2bÞnþ1 coshðbxÞ n odd 0 1 ! n n1 2 2 X 1 X 1 2s 2sþ1 x an;2s coshðbxÞ þ @ x an;2sþ1 A sinhðbxÞ A ¼ 4b s¼0 4b s¼0 sinhðbxÞ n odd n! nþ1 coshðbxÞ n even 2ð2bÞ
Substituting back into Eq. (18) yields the following result: I nss ¼ Y ns sinhðbnÞ þ Z ns coshðbnÞ; I ncs ¼ Y nc sinhðbnÞ þ Z nc coshðbnÞ I sc ¼ K ns sinhðbnÞ þ P ns coshðbnÞ; I cc ¼ K nc sinhðbnÞ þ P nc coshðbnÞ
where
3948
G.E. Cossali / International Journal of Heat and Mass Transfer 50 (2007) 3943–3948
! n 2 1 X 1 þ ð1Þn n! n2s ¼ x an;n2s þ ; nþ1 4b s¼0 2 2ð2bÞ 0 1 n1 n 2 nþ1 X 1 x 1 ð1Þ n! n n2sþ1 A x an;n2sþ1 þ Z s ¼ @ 2ðn þ 1Þ 2 4b s¼0 2ð2bÞnþ1 0 1 n1 n 2 nþ1 X 1 x 1 ð1Þ n! A; xn2sþ1 an;n2sþ1 Y nc ¼ @ nþ1 2ðn þ 1Þ 2 4b s¼0 2ð2bÞ ! n n 2 1 X 1 þ ð1Þ n! n n2s Zc ¼ x an;n2s 4b s¼0 2 2ð2bÞnþ1 0 1 n1 n 2 nþ1 X 1 x 1 ð1Þ n! A; K ns ¼ @ xn2sþ1 an;n2sþ1 þ nþ1 2ðn þ 1Þ 2 4b s¼0 2ð2bÞ ! n n 2 1 X 1 þ ð1Þ n! P ns ¼ xn2s an;n2s 2 4b s¼0 2ð2bÞnþ1 ! n n 2 1 X 1 þ ð1Þ n! n n2s Kc ¼ x an;n2s þ ; 2 4b s¼0 2ð2bÞnþ1 0 1 n1 n 2 nþ1 X 1 x n! 1 ð1Þ A P nc ¼ @ xn2sþ1 an;n2sþ1 þ þ 4b s¼0 2ðn þ 1Þ 2ð2bÞnþ1 2 Y ns
n
ð1Þ n! and an;k ¼ k!ð2bÞ nk .
References [1] A.I. Nakorchevkii, Dynamics of ground heat storage and choice of rational solutions, J. Eng. Phys. Thermophys. 77 (4) (2004) 688–699. [2] K.A. Antonopoulos, F. Democritou, Transient conduction under cosine temperature perturbations, Appl. Energy (50) (1995) 233–246. [3] S.A. Adamovsky, A.A. Minakov, C. Schick, Scanning microcalorimetry at high cooling rate, Thermochim. Acta 403 (2003) 55–63. [4] P.T. Ireland, T.V. Jones, Detailed measurements of heat transfer on and around a pedestal in fully developed passage flow, in: Proceedings of the 8th International Heat Transfer Conference Proceedings, vol. 3, 1986, pp. 975–980. [5] R.J. Butler, J.W. Baughn, The effect of the thermal boundary conditions on transient method heat transfer measurements on a. at plate with a laminar boundary layer, J. Heat Transfer 118 (1996) 831– 837.
[6] F. De Monte, Transient heat conduction in one-dimensional composite slab. A ‘natural’ analytic approach, Int. J. Heat Mass Transfer 43 (2000) 3607–3619. [7] B-L. Wang, Z-H. Tian, Application of finite element-finite difference method to the determination of transient temperature field in functionally graded materials, Finite Element Anal. Design 41 (2005) 335–349. [8] A. Sutradhar, G.H. Paulino, L.J. Gray, Transient heat conduction in homogeneous and non-homogeneous materials by the Laplace transform Galerkin boundary element method, Eng. Anal. Boundary Elements 26 (2002) 119–132. [9] H. Asan, Y.S. Sancaktar, Effects of wall’s thermophysical properties on time lag and decrement factor, Energy Build. 28 (1998) 159–166. [10] H. Asan, Numerical computation of time lags and decremental factors for different building materials, Build. Environ. 41 (2006) 615– 620. [11] A.A. Minakov, S.A. Adamovsky, C. Schick, Advanced two-channel a.c. calorimeter for simultaneous measurements of complex heat capacity and complex thermal conductivity, Thermochim. Acta 403 (2003) 89–103. [12] U.G. Jonsson, O. Andersson, Investigation of the low- and highfrequency response of 3-x sensors used in dynamic heat capacity measurements, Meas. Sci. Technol. 9 (1998) 1873–1885. [13] G.E. Cossali, A Fourier transform based data reduction method for the evaluation of the local convective heat transfer coefficient, Int. J. Heat Mass Transfer 47 (2004) 21–30. [14] Y. Chen, J. Zhou, J.D. Spitler, Verification for transient heat conduction calculation of multilayer building constructions, Energy Build. 38 (2006) 340–348. [15] F. De Monte, An analytic approach to the unsteady heat conduction processes in one-dimensional composite media, Int. J. Heat Mass Transfer 45 (2002) 1333–1343. [16] F. Alhama, J.F. Lopez-Sanchez, C.F. Gonzales-Fernandez, Heat conduction through a multilayered wall with variable boundary conditions, Energy 22 (8) (1997) 797–803. [17] C-H. Chen, J.-S. Chiou, Periodic heat transfer in a vertical plate fin cooled by a forced convective flow, Int. J. Heat Mass Transfer 39 (1996) 429–435. [18] Y. Chen, S. Wang, A new procedure for calculating periodic response factors based on frequency domain regression method, Int. J. Thermal Sci. 44 (2005) 382–392. [19] S. Zubair, M. Aslam Chaudhry, Heat conduction in a semi-infinite solid subject to steady and non-steady periodic-type surface heat fluxes, Int. J. Heat Mass Transfer 38 (1995) 3393–3399. [20] D. Maillet, S. Andre`, J.C. Batsale, A. Degiovanni, C. Moyne, Thermal Quadrupoles, Wiley, 2000.
Dynamic response of a non-homogeneous 1D slab under periodic thermal excitation G.E. Cossali * Facolta` di Ingegneria, Universita` di Bergamo, via Marconi 6, 24044 Dalmine (BG), Italy Received 6 June 2006 Available online 30 March 2007
Abstract The periodic heat conduction in a non-homogeneous slab is studied and an approximate analytic solution, valid for weakly non homogeneous slabs, is found. An exact solution for a well defined non-homogeneity is also found and compared to the approximate solution to estimate the accuracy of the result. The solution allows to extend the method of thermal quadrupoles to weakly non-homogeneous 1D-slabs. Ó 2007 Elsevier Ltd. All rights reserved.
1. Introduction Transient heat conduction in homogeneous materials is a deeply analysed topic, also for the vast consequences in many applied fields, like energy storage [1], laser heating, thermal behaviour of building materials [2], but also for measurement techniques like microcalorimetry [3] or transient techniques for measuring convective heat transfer coefficient [4,5]. Many different techniques were developed to solve transient one-dimensional problems also for the more complex case of composite walls [6]. The case of transient heating of a non-homogeneous materials is more complex and numerical approaches were used to find solutions for a certain class of materials [7,8]. Periodic heat conduction is a relatively less investigated topic, although the implications for many applied fields are significant, for example for building walls and materials [9,10] or again for techniques to measure conductivity and heat capacity like a.c. calorimetry [11] or 3-x methods [12] and for local convective heat transfer coefficient measurements [13]. The periodic conduction in a 1D homogeneous slab is a relatively simple problem but recently many studies, dealing
with it, can be found in the literature [16–19], particularly related to the hyperbolic version of the heat equation, as the wave like behaviour can be better evidenced in a simple geometry. However, the extension of the transient conduction analysis to non-homogeneous material is not straightforward, and fewer works can be found on this subject [14,15], whereas periodic conduction in non-homogeneous material does not seems to have been sufficiently considered, although non-homogeneous material are becoming everyday more widespread in many engineering applications. The paper presents a solution of the periodic conduction in a non-homogeneous 1D slab under the condition of weak non-homogeneity (i.e. limited variation of the material properties with position) and an assessment of the accuracy of the method by comparison with an exact solution. The solution allows a relatively simple extension of the thermal quadrupole method to this class of materials. 2. Basic equations Consider a 1D slab made of a non-homogeneous material, the energy equation can then be written as qc
*
Tel.: +39 03520 52309; fax: +39 03556 2779. E-mail address: [email protected]
0017-9310/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2007.01.038
oT oq ¼ ot ox
and making use of Fourier law: q ¼ k oT we obtain: ox
ð1Þ
3944
d2
G.E. Cossali / International Journal of Heat and Mass Transfer 50 (2007) 3943–3948
oT oT o2 T ¼ cðnÞ þ aðnÞ 2 ot on on
ð2Þ
ok ; and aðnÞ ¼ qck . The temperature where n ¼ dx ; cðnÞ ¼ qc1 on and heat flux fields can be split into time average (T a ; qa ) and fluctuating ( Te , ~ q) components. Introducing the Fourier transform of the latter: Z þ1 T ðx; tÞ ¼ T a ðxÞ þ Sðx; xÞeixt dx; 1 Z þ1 Qðx; xÞeixt dx qðx; tÞ ¼ qa ðxÞ þ 1
Eq. (2) becomes: ixd2 S ¼ cðnÞS 0 þ aðnÞS 00
ð3Þ
0 ¼ cðnÞT 0a þ aðnÞT 00a
ð4Þ
where apex means derivation respect to n. The case of homogeneous slab is found setting c ¼ 0 and a ¼ a0 ¼ const: and Eq. (3) becomes: S 00H b2 S H ¼ 0 qffiffiffiffiffiffiffi 2 with b ¼ ixd , the solution is easily found under the more a0 convenient form: S H ¼ S H;0 cosh ðbnÞ þ
S 0H;0 sinh ðbnÞ b
ð5Þ
with S H;0 ¼ S H ð0Þ and S 0H;0 ¼ S 0H ð0Þ. The general solution of Eq. (4) is easily found by double integration: Z n R f cðvÞ dv Ta ¼ A þ B e 0 aðvÞ df: 0
mðnÞ ¼ S H0 m1 ðnÞ þ
S 0H0 m2 ðnÞ b
ð7Þ
cðnÞ aðnÞ 2 m1 ðnÞ ¼ b sinh ðbnÞ þ b cosh ðbnÞ ; a0 a0 cðnÞ aðnÞ 2 m2 ðnÞ ¼ b coshðbnÞ þ b sinhðbnÞ a0 a0 are known functions. The solution of (6) is then: Z 1 1 n s ¼ s0 coshðbnÞ þ s00 sinhðbnÞ mðzÞ sinh½bðz nÞdz b b 0 and choosing sð0Þ ¼ 0, s0 ð0Þ ¼ 0, we obtain: Z 1 n S ¼ S 0 cosh ðbnÞ m1 ðzÞ sinh½bðz nÞdz b 0 Z S0 1 n þ 0 sinh ðbnÞ m2 ðzÞ sinh½bðz nÞdz b 0 b
ð8Þ
with S 0 ¼ Sð0Þ and S 00 ¼ S 0 ð0Þ. The linearisation that led to Eq. (6) introduces approximation whose entity depends on the magnitude of the slab thermal properties variations. To assess the validity of the approximation a comparison to an exact solution of Eq. (3), for a particular distribution of the slab properties, is needed. 4. An exact solution Consider the following particular form for the dependence on position of the material properties: qc ¼ const; k ¼ Kðn0 nÞ
2
ð9Þ
Eq. (3) becomes qcxd2 S¼0 K and after the change of qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi variable: g ¼ ðn n0 Þ, and defin4xqcd2 1 ing k0 ¼ 2 1 þ i K , an Euler equation is promptly recovered: 1 2 0 2 00 k0 S ¼ 0 g S þ 2gS þ ð10Þ 4 2
ðn0 nÞ S 00 þ 2ðn n0 ÞS 0 i
3. The weakly non-homogeneous slab Consider a material characterised by variations of its properties with position that are small compared to the average value, i.e.: k ¼ k 0 þ jðnÞ;
jjðnÞj jk 0 j
with
0
0
jj ðnÞj ¼ jk ðnÞj jkðnÞj a ¼ a0 þ aðnÞ;
with
jaðnÞj ja0 j 0
from which it follows: j ac j ¼ j kk j 1. Suppose now that the solution of Eq. (3) can be written as a perturbation of the solution SH valid for homogeneous materials (Eq. (5)), i.e. S ¼ SH þ s
with
s ¼ oðS H Þ
then substituting into Eq. (3) and preserving only the terms up to the first order: s00
ixd 2 cðnÞ 0 aðnÞ 00 s¼ S S ¼ mðnÞ a0 H a0 H a0
where
ð6Þ
where now the apex means derivation respect to g. The solution of this equation can be written in the convenient form: 1 1 n0 n 2 1 S¼ S 0 sinhðzÞ þ k0 coshðzÞ S 00 n0 sinhðzÞ k0 n0 2 ð11Þ with
n n z ¼ k0 ln 0 n0
and again: Sð0Þ ¼ S 0 , S 0 ð0Þ ¼ S 00 . The exact analytic solution can now be used to assess the validity of the approximate solution (8).
G.E. Cossali / International Journal of Heat and Mass Transfer 50 (2007) 3943–3948
(
5. Assessment of the perturbative solution F 1;ap Consider again the material described by (9) with the further condition n0 1, then, after defining: k 0 ¼ Kn20 we obtain:
F 2;ap
) ! N 1 ðnÞ M 1 ðnÞ ¼ 1þ 2 sinh ðbnÞ cosh ðbnÞ þ n0 n20 ( ) " # 1 M 2 ðnÞ N 2 ðnÞ sinh ðbnÞ 1 þ ¼ þ 2 coshðbnÞ b n20 n0
jðnÞ ¼ k k 0 ¼ Kn2 2Kn0 n k 0 ðxÞ ¼ 2Kðn n0 Þ
ð13Þ with the functions:
and the condition: jk 0 j jkj yields: 2 jðn n0 Þj that is obviously verified as 0 6 n 6 1 and n0 1. Moreover a0 ¼
k 0 Kn20 ¼ ; qc qc
cðnÞ 2ðn n0 Þ ¼ ; a0 n20
3945
aðnÞ n2 2n0 n ¼ a0 n20
The weakly non-homogeneous approximation holds and the functions m1 and m2 (see Eq. (7)) now become: ( ) 2ðn n0 Þ n2 2n0 n 2 m1 ðnÞ ¼ b sinh ðbnÞ þ b coshðbnÞ ; n20 n20 ( ) 2ðn n0 Þ n2 2n0 n 2 b cosh ðbnÞ þ b sinhðbnÞ m2 ðnÞ ¼ n20 n20 qffiffiffiffiffiffiffiffi 2 with b ¼ i xd and the following integrals: a0 Z n m1 ðzÞ sinh ðb½z xÞ dz I 1 ðnÞ ¼
1 1 n0 n 2 1 n0 n F 1;ex ¼ sinh k0 ln k0 n0 2 n0 n0 n þk0 cosh k0 ln ð14Þ n0 1 n0 n 2 n0 n0 n F 2;ex ¼ sinh k0 ln n0 n0 k0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi respectively, with now k0 ¼ 12 1 þ 4b2 n20 . Fig. 1 reports a comparison for a fixed value of x, having introduced the non-homogeneity magnitude e defined as e¼
k max k min jð1 2n0 Þj ¼ k0 n20
0
I 2 ðnÞ ¼
Z
n
m2 ðzÞ sinh ðb½z xÞ dz
0
must be evaluated to find the solution in an explicit form. In Appendix, the integrals are calculated, finding: b fM 1 sinhðbnÞ þ N 1 cosh ðbnÞg n20 b I 2 ðnÞ ¼ 2 fM 2 sinhðbnÞ þ N 2 cosh ðbnÞg n0 I 1 ðnÞ ¼
with n0 1 n b nn 1 þ n 0 bn2 n3 ; N 1 ¼ 0 n2 6 2 4 2b 4b 2 1 n0 1 2 1 n0 b 2 b 3 n n n M2 ¼ 2 n n ; N2 ¼ 4 4b 2 6 2 4b M1 ¼
The approximate solution then becomes: ( ) ! N 1 ðnÞ M 1 ðnÞ S ¼ S H0 1þ sinhðbnÞ cosh ðbnÞ þ n20 n20 ( ) " # S 0H0 M 2 ðnÞ N 2 ðnÞ sinh ðbnÞ 1 þ þ þ 2 cosh ðbnÞ b n20 n0 ð12Þ The exact (11) and approximate (12) solutions can now be compared by comparing the functions:
Fig. 1. Relative discrepancies between the approximated and the exact values of: (a) the function F1 and (b) the function F2 versus the position along the slab, for x ¼ 1.
3946
G.E. Cossali / International Journal of Heat and Mass Transfer 50 (2007) 3943–3948
It appears that the maximum discrepancy occurs at n ¼ 1, thus it is interesting to compare the two solutions (by comparing the functions (13) and (14)) at the fixed position n ¼ 1 for different values of the exciting frequency. To this end, observe that by introducing the parameter ~e ¼ n1 0 ¼ e þ oðeÞ, the coefficients of the hyperbolic functions in F j;a evaluated for n ¼ 1 can be expressed as polynomials of second order in ~e: ~e ~e2 N1 ¼ 1 2 4 n20 ~e 1 ~e2 1 b M1 b þ C 12;a ¼ 2 ¼ 2 2b 3 2 b n0 2 ~e ~e 1 N2 1 1 C 21;a ¼ ¼ þ b n20 2 2 2b2 3 ! ~e ~e2 1 1 M2 1 1þ 2 ¼ 1 C 22;a ¼ þ1 b b 2 4 b2 n0
Moreover, defining: u ¼ k0 ln n0n1 b the functions F j;e 0 can be written (for n ¼ 1) as 1 1 sinh ðuÞ þ cosh ðuÞ cosh ðbÞ F 1;ex ¼ ð1 ~eÞ 2 2k0 1 1 2 þð1 ~eÞ coshðuÞ þ sinh ðuÞ sinh ðbÞ 2k0 C 11;a ¼ 1 þ
12
F 2;ex
ð1 ~eÞ ¼ ~ek0
sinh ðuÞ cosh ðbÞ 1
ð1 ~eÞ 2 þ cosh ðuÞ sinh ðbÞ ~ek0
!
Fig. 2. Relative discrepancies between the approximated and the exact values of: (a) the function F1 and (b) the function F2 versus nondimensional frequency (b) for n ¼ 1.
and the corresponding coefficients of the hyperbolic functions are (up to Second order in ~e): ~e ~e2 C 11;e ¼ 1 þ ð1 þ b2 Þ þ oð~e2 Þ 8 2 ~e 1 ~e2 7b 1 b þ C 12;e ¼ þ oð~e2 Þ 2 b 4 3 2b ~e ~e2 7b2 1 C 21;e ¼ þ þ þ oð~e2 Þ 2 4 3b2 2b2 ! ~e 1 2 b4 þ 3b2 1 1 2 1 þ ~e þ oð~e Þ C 22;e ¼ b 2 8 b2
ð15Þ 6. The thermal quadrupole approach for the non-homogeneous slab
a comparison with those relative to the approximate solution shows that the error is of order ~e2 , i.e. of order e2 ~e2 ð3 þ b2 Þ þ oð~e2 Þ 8 3~e2 1 b þ oð~e2 Þ ¼ 4 2b ~e2 1 ¼ 3 þ oð~e2 Þ 4 2b2 ~e2 b2 þ 4 ¼ þ oð~e2 Þ 2b 4
DC 11 ¼ C 11;a C 11;e ¼ DC 12 ¼ C 12;a C 12;e DC 21 ¼ C 21;a C 21;e DC 22 ¼ C 22;a C 22;e
Fig. 2 shows the comparison of the two solutions, the discrepancy between the approximate and the exact one remains limited also for relatively large (0.2) values of e. However, relative errors (DC jk =C jk;e ) are expected to increase for very large and very small values of the frequency x, as it stems from Eqs. (15) and (16).
ð16Þ
The thermal quadrupole approach [20] gives a compact and elegant way to solve transient problems for 1D homogeneous slabs and the application to multilayered slabs is straightforward. The previous results allow to extend the thermal quadrupole approach also to 1D non-homogeneous materials. Defining S Wðx; xÞ ¼ Q the following relation, that is the basis of the thermal quadrupole method, holds for a homogeneous slab: Wðd; xÞ ¼ Mh Wð0; xÞ where
G.E. Cossali / International Journal of Heat and Mass Transfer 50 (2007) 3943–3948
" Mh ¼
cosh ðbnÞ
1 sinh ðbnÞ kb
kb sinh ðbnÞ
cosh ðbnÞ
#
and using partial integration, the following recurrence relation is found:
For the non-homogeneous slab a similar equation can now be written, in fact: Z 1 n S ¼ S 0 cosh ðbnÞ m1 ðzÞ sinh½bðz nÞ dz b 0 Z Q0 1 n sinh ðbnÞ m2 ðzÞ sinh½bðz nÞ dz b 0 kb Z n Q ¼ S 0 kb sinh ðbnÞ k m1 ðzÞ cosh½bðz nÞ dz þ Q0 cosh ðbnÞ þ
1 b
Z
0 n
m2 ðzÞ cosh½bðz nÞ dz
"
coshðbnÞ
I n ¼
e2bx n n I x 2b n1 2b
where: Z I ¼ 0
x
expð2bzÞ dz ¼
0
e2bx 1 2b
and the closed form: I n ¼
nk n e2bx X ð1Þ n! n! nþ1 xk þ ð1Þ nþ1 2b k¼0 k!ð2bÞnk ð2bÞ
ð17Þ
0
and then defining: Mnh ¼
3947
kb1 sinhðbnÞ
#
kb sinhðbnÞ coshðbnÞ 2 Rn b1 0 m1 ðzÞ sinh½bðz nÞdz þ4 Rn k 0 m1 ðzÞ cosh½bðz nÞdz " 1 # b I 1 kb12 I 2 ¼ Mh þ kI 3 b1 I 4
1 kb2 1 b
Rn
Rn 0
0
m2 ðzÞ sinh½bðz nÞdz
m2 ðzÞ cosh½bðz nÞdz
3 5
is promptly found. Consider now the four integrals: Z x n I ss ¼ zn sinh ðbzÞ sinh ðbðz xÞÞ dz Z0 x n I cs ¼ zn cosh ðbzÞ sinh ðbðz xÞÞ dz 0 Z x I sc ¼ zn sinh ðbzÞ cosh ðbðz xÞÞ dz 0 Z x I cc ¼ zn cosh ðbzÞ cosh ðbðz xÞÞ dz 0
the following equation is obtained: Wðd; xÞ ¼ Mnh Wð0; xÞ
a simple calculation gives:
If the functions cðnÞ and aðnÞ can be expressed as polynomial in n, then the integrals Ij can be written in a closed form (see Appendix). These results allow to extend the use of the thermal quadrupole method to 1-d weakly non-homogeneous slabs.
xnþ1 xnþ1 coshðbxÞ; I ncs ¼ A sinhðbxÞ; 2ðn þ 1Þ 2ðn þ 1Þ xnþ1 xnþ1 I nsc ¼ A þ sinhðbxÞ; I ncc ¼ Aþ þ coshðbxÞ 2ðn þ 1Þ 2ðn þ 1Þ
7. Conclusions The thermal response of a weakly non-homogeneous slab to periodic thermal excitation was analysed through a Fourier transform approach and an analytic solution of the transformed problem, valid in the limit of small nonhomogeneity, was found. The solution can be written in closed form whenever the dependence of the wall thermal characteristics on the position can be represented in polynomial form. To assess the validity of the approximation, a comparison with an exact solution of the transformed problem is reported, showing that accuracy is of second order in the non-homogeneity parameter e. The solution allows also to extend the method of thermal quadrupoles to the weakly non-homogeneous slab thus extending the capability of the method. Appendix Defining the integrals: Z x In ¼ zn expð2bzÞ dz 0
I nss ¼ Aþ
ð18Þ where 1 bx A ¼ ½ebx I þ n e In 4 and using Eq. (17), after few manipulations, yields: 0 1 ! n n1 2 2 X 1 X 1 Aþ ¼ xn2s an;n2s sinhðbxÞ þ @ xn2sþ1 an;n2sþ1 A coshðbxÞ 4b s¼0 4b s¼0 sinhðbxÞ n even n! þ 2ð2bÞnþ1 coshðbxÞ n odd 0 1 ! n n1 2 2 X 1 X 1 2s 2sþ1 x an;2s coshðbxÞ þ @ x an;2sþ1 A sinhðbxÞ A ¼ 4b s¼0 4b s¼0 sinhðbxÞ n odd n! nþ1 coshðbxÞ n even 2ð2bÞ
Substituting back into Eq. (18) yields the following result: I nss ¼ Y ns sinhðbnÞ þ Z ns coshðbnÞ; I ncs ¼ Y nc sinhðbnÞ þ Z nc coshðbnÞ I sc ¼ K ns sinhðbnÞ þ P ns coshðbnÞ; I cc ¼ K nc sinhðbnÞ þ P nc coshðbnÞ
where
3948
G.E. Cossali / International Journal of Heat and Mass Transfer 50 (2007) 3943–3948
! n 2 1 X 1 þ ð1Þn n! n2s ¼ x an;n2s þ ; nþ1 4b s¼0 2 2ð2bÞ 0 1 n1 n 2 nþ1 X 1 x 1 ð1Þ n! n n2sþ1 A x an;n2sþ1 þ Z s ¼ @ 2ðn þ 1Þ 2 4b s¼0 2ð2bÞnþ1 0 1 n1 n 2 nþ1 X 1 x 1 ð1Þ n! A; xn2sþ1 an;n2sþ1 Y nc ¼ @ nþ1 2ðn þ 1Þ 2 4b s¼0 2ð2bÞ ! n n 2 1 X 1 þ ð1Þ n! n n2s Zc ¼ x an;n2s 4b s¼0 2 2ð2bÞnþ1 0 1 n1 n 2 nþ1 X 1 x 1 ð1Þ n! A; K ns ¼ @ xn2sþ1 an;n2sþ1 þ nþ1 2ðn þ 1Þ 2 4b s¼0 2ð2bÞ ! n n 2 1 X 1 þ ð1Þ n! P ns ¼ xn2s an;n2s 2 4b s¼0 2ð2bÞnþ1 ! n n 2 1 X 1 þ ð1Þ n! n n2s Kc ¼ x an;n2s þ ; 2 4b s¼0 2ð2bÞnþ1 0 1 n1 n 2 nþ1 X 1 x n! 1 ð1Þ A P nc ¼ @ xn2sþ1 an;n2sþ1 þ þ 4b s¼0 2ðn þ 1Þ 2ð2bÞnþ1 2 Y ns
n
ð1Þ n! and an;k ¼ k!ð2bÞ nk .
References [1] A.I. Nakorchevkii, Dynamics of ground heat storage and choice of rational solutions, J. Eng. Phys. Thermophys. 77 (4) (2004) 688–699. [2] K.A. Antonopoulos, F. Democritou, Transient conduction under cosine temperature perturbations, Appl. Energy (50) (1995) 233–246. [3] S.A. Adamovsky, A.A. Minakov, C. Schick, Scanning microcalorimetry at high cooling rate, Thermochim. Acta 403 (2003) 55–63. [4] P.T. Ireland, T.V. Jones, Detailed measurements of heat transfer on and around a pedestal in fully developed passage flow, in: Proceedings of the 8th International Heat Transfer Conference Proceedings, vol. 3, 1986, pp. 975–980. [5] R.J. Butler, J.W. Baughn, The effect of the thermal boundary conditions on transient method heat transfer measurements on a. at plate with a laminar boundary layer, J. Heat Transfer 118 (1996) 831– 837.
[6] F. De Monte, Transient heat conduction in one-dimensional composite slab. A ‘natural’ analytic approach, Int. J. Heat Mass Transfer 43 (2000) 3607–3619. [7] B-L. Wang, Z-H. Tian, Application of finite element-finite difference method to the determination of transient temperature field in functionally graded materials, Finite Element Anal. Design 41 (2005) 335–349. [8] A. Sutradhar, G.H. Paulino, L.J. Gray, Transient heat conduction in homogeneous and non-homogeneous materials by the Laplace transform Galerkin boundary element method, Eng. Anal. Boundary Elements 26 (2002) 119–132. [9] H. Asan, Y.S. Sancaktar, Effects of wall’s thermophysical properties on time lag and decrement factor, Energy Build. 28 (1998) 159–166. [10] H. Asan, Numerical computation of time lags and decremental factors for different building materials, Build. Environ. 41 (2006) 615– 620. [11] A.A. Minakov, S.A. Adamovsky, C. Schick, Advanced two-channel a.c. calorimeter for simultaneous measurements of complex heat capacity and complex thermal conductivity, Thermochim. Acta 403 (2003) 89–103. [12] U.G. Jonsson, O. Andersson, Investigation of the low- and highfrequency response of 3-x sensors used in dynamic heat capacity measurements, Meas. Sci. Technol. 9 (1998) 1873–1885. [13] G.E. Cossali, A Fourier transform based data reduction method for the evaluation of the local convective heat transfer coefficient, Int. J. Heat Mass Transfer 47 (2004) 21–30. [14] Y. Chen, J. Zhou, J.D. Spitler, Verification for transient heat conduction calculation of multilayer building constructions, Energy Build. 38 (2006) 340–348. [15] F. De Monte, An analytic approach to the unsteady heat conduction processes in one-dimensional composite media, Int. J. Heat Mass Transfer 45 (2002) 1333–1343. [16] F. Alhama, J.F. Lopez-Sanchez, C.F. Gonzales-Fernandez, Heat conduction through a multilayered wall with variable boundary conditions, Energy 22 (8) (1997) 797–803. [17] C-H. Chen, J.-S. Chiou, Periodic heat transfer in a vertical plate fin cooled by a forced convective flow, Int. J. Heat Mass Transfer 39 (1996) 429–435. [18] Y. Chen, S. Wang, A new procedure for calculating periodic response factors based on frequency domain regression method, Int. J. Thermal Sci. 44 (2005) 382–392. [19] S. Zubair, M. Aslam Chaudhry, Heat conduction in a semi-infinite solid subject to steady and non-steady periodic-type surface heat fluxes, Int. J. Heat Mass Transfer 38 (1995) 3393–3399. [20] D. Maillet, S. Andre`, J.C. Batsale, A. Degiovanni, C. Moyne, Thermal Quadrupoles, Wiley, 2000.