A general modal-based numerical simulation of transient heat conduction in a one-dimensional homogeneous slab

Energy and Buildings 42 (2010) 2309–2322

Contents lists available at ScienceDirect

Energy and Buildings journal homepage: www.elsevier.com/locate/enbuild

A general modal-based numerical simulation of transient heat conduction in a one-dimensional homogeneous slab G. Lefebvre Centre d’Études et de Recherches en Thermique, Environnement et Systèmes CERTES EA3481, Université Paris-Est, 94010 Créteil, France

a r t i c l e

i n f o

Article history: Received 1 March 2010 Received in revised form 21 June 2010 Accepted 20 July 2010 Keywords: Heat conduction Modal Eigen function Eigen value Orthogonal expansion Proper order decomposition Model size reduction Diffusion process

a b s t r a c t We present an analytical general solution of the well-known heat conduction problem in a onedimensional homogeneous slab by using an expansion of the only delayed part of the temperature field function on an infinite eigen functions set. We propose to split the searched solution of the general problem in an instantaneous part and a delayed one which, with whatever boundary conditions combination, verifies Dirichlet conditions and then avoids the Gibbs problem. A space sampling and a reduction of the eigen basis produce a useful and usable numerical model which may be as accurate as required. The results we present may be used directly to numerically calculate solutions to a given problem and may be reused as a basis for more complex problems as multilayer walls or whole thermal systems studies. © 2010 Elsevier B.V. All rights reserved.

1. Introduction

2. Hypothesis and governing equations

Analyzing and quantifying heat conduction through a plain homogeneous thin slab seems simple nowadays, thanks to the huge number of published analytical (method of the separation of variables, symbolic transforms, Green functions, etc.) and numerical solutions (finite differences or elements, control volume method, coefficients of transfer functions, etc.) (see among many others [1–6]). This paper presents a simple and, in our mind, pertinent analytical solution of the well-known conductive heat transfer problem in a one-dimensional homogeneous slab. However the obtained solution is very general as it is valid for any combination of boundary, heat source and initial conditions. Published works use the spectral characteristics of the heat transfer and boundary condition equations to solve the heat transfer problem [7,8]. The proposed way uses an expansion of the purely dynamic part of the temperature field function on an infinite set of eigen functions which constitutes a basis of the functional space of possible solutions. A space sampling and a reduction of the eigen basis produce a useful and usable numerical model which may be as accurate as required.

We consider a solid homogeneous cuboid (a slab) with one dimension L small (the thickness) versus the two others (Fig. 1). The expected temperature range allows us to suppose that the solid and boundary characteristics (heat capacity c, heat conductivity k, thermal diffusivity a and the superficial heat transfer coefficients hi ) are constant. The thermal state described by the temperature field, only depends on the time t and the space variable x in the thickness direction. The boundary condition which may be applied independently on sides a and b of the slab is the sum of a convective term (Newton–Robin relation including a heat transfer coefficient hi ) and an absorbed heat flux term ϕi , each being positive when heat enters the slab. The proposed expression of the right-hand side of the equations (a 3.5 type following the typology proposed in [9,10]) may represent all classical boundary condition types, by choosing particular values for the adequate parameters and variables:

Abbreviations: NSM, not split modal model; SM, split modal model. E-mail address: [email protected]. 0378-7788/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.enbuild.2010.07.024

- 1st type: Dirichlet condition: hi = ∞; ϕi = 0. - 2nd type: Neumann condition: hi = 0. - 3rd type: Newton–Robin condition: ϕi = 0. The mathematical model which describes the heat conduction problem is a set of four Eq. (1) the heat equation, (2) and (3) two boundary condition equations, and (4) the initial condition equation. Eq. (1) is a partial differential equation with a first order

2310

G. Lefebvre / Energy and Buildings 42 (2010) 2309–2322

Nomenclature a A(t) ai , bi B(t)

thermal diffusivity (m2 s−1 ) coefficient of the sliding temperature T◦ expression (K m−1 ) parameters of eigen function i coefficient of the sliding temperature T◦ expression (K) Biot numbers on side a and b (.) coefficient of the state space equation specific heat capacitance (J kg−1 K−1 ) integration constants (T¯ expression) (K) constants used to express T¯

Ba , Bb Bij c c1 , c2 c1. , c2. D(t), E(t) equivalent temperature gradient (boundary conditions) (K m−1 ) f(x) initial temperature field (K) G(x,t) temperature field (K) ha , hb convective heat exchange coefficients on sides a and b (W m−2 K−1 ) k heat conductivity of the material (W m−1 K−1 ) L slab thickness (m) L heat operator (.) Ni norm of an eigen function P Biot numbers product (.) p(x,t) volumic heat source (W m−3 ) pk elementary space function of the heat source expansion (W m−3 ) R reduction order S Biot numbers sum (.) Sk (x) elementary space function of the expanded sliding temperature T¯ S t t T(x,t) Tˆ T¯

formal solution operator (K) time (s) time step (s) space and time temperature function (K) complement of T¯ in the general solution (K) temperature field when there is no heat capacitance (K) T◦ temperature field when there is neither heat capacitance nor heat source (K) T+ complement of T◦ in the general solution (K) T˜ Approximate expression of the general solution (K) initial time (s) t0 T functional space of functions defined on]0,L[(K) Uk (t) elementary time function of the expanded sliding temperature = solicitation elementary time function of the equivalent heat uk (t) source expansion (.) Vi (x) eigen function X Spatial coordinate in the slab thickness direction (m) Xi (t), xi (t) modal states Z(t) source term expressed as a temperature variability (K s−1 ) Z* (t) fictitious source term in the T+ equation (K s−1 ) zk (x) elementary space function of the equivalent heat source expansion (K s−1 )

 i ,

material density (kg m−3 ) time constant of eigen element i (s) integration variables

Indices 0 a b p R

initial left side right side steady state reduction order

Superscripts ∞ steady state * instantaneous (temperature field) + Delayed (temperature field)

differential versus time, requiring one equation (4), and a second order differential versus space, requiring two Eqs. (2) and (3) to determine the unknown parameters which appear when integrating. The heat equation is classically obtained by using the Fourier conduction equation in a heat balance equation derived from the thermodynamics first principal. Eventual heat sources (p > 0) or sinks (p < 0) are described by the right-hand side term p in Eq. (1). ∂T ∂2 T (x, t) = a 2 (x, t) + Z(x, t) ∀x ∈ ]0, L[ ∂t ∂x

∀t ≥ t0

(1)

−

∂T Ba T (x, t) = D(t) (x, t) + L ∂x

x = 0 ∀t ≥ t0

(2)

+

∂T B (x, t) + b T (x, t) = E(t) L ∂x

x=L

∀t ≥ t0

(3)

T (x, t) = f (x) ∀x ∈ ]0, L[

t = t0

(4)

where Z(x, t) = p(x, t)/c is a source of temperature variability, D(t) = ϕa /k + Ba Ta /L and E(t) = ϕb /k + Bb Tb /L are equivalent boundary temperature gradient. Ta , Tb , ϕa , ϕb and p are solicitations which thermally stimulate the slab and determine, with the initial condition, its future thermal state. The objective is to find an expression of the temperature T(x,t) in function of x, t, the parameters c, k, ha , hb and the five solicitations. Ba = ha L/k and Bb = hb L/k are the dimensionless Biot numbers on respective sides a and b and a = k/c is the thermal diffusivity.

Greek letters ˛i parameter of eigen function i ε temperature error (K) ϕa (t), ϕb (t) surfacic heat flux absorbed on sides a and b (W m−2 ) i Eigen value (s−1 ) Fig. 1. Homogeneous 1D slab with boundary conditions.

G. Lefebvre / Energy and Buildings 42 (2010) 2309–2322

2311

3. Solving with the modal method: principals

3.2. Norm of an eigen function

The mathematical fundamentals on which relies the solution are:



- a space function of x ∈ ]0, L[ as T(x, t) may be seen as an element of a functional space T. - it is assumed that an infinite denumerable set of linearly independent functions of T is a basis of T. - a basis of this functional space will be used to express at any time a part of the temperature field as a linear expansion of the basis functions. - the part which will be expanded is chosen to present fast numerical convergence of the resulting truncated series. Among all the possibilities and in using the same idea as Fourier mentioned in 1822 [11], the basis functions are chosen in order to have a physical meaning and to produce a pertinent model. The basis is made of the eigen functions of the heat operator L{f}= a ∂ 2 f/∂ x2 which verify boundary conditions (2) and (3). If i is an eigen value of L and Vi is the corresponding eigen function, we obtain: L{Vi } = a

∂2 Vi = i Vi ∂x2

(5)

In our particular case of a self-adjoint operator, all associated eigen values are real, distinct and simple [12]. We know that the solution of this equation, called the “Helmholtz equation”, is: Vi (x) = ai cos

˛ x i

L

+ bi sin

˛ x i

L

∀x ∈ [0, L]

(6)

where ai , bi and ˛i are arbitrary constants which will be determined later thanks to initial and boundary conditions applied on the slab. The ith eigen value i is related to the constants value:

 ˛ 2

i = −a

i

f1 , f2



1 = L

L

f1 (x) f2 (x) dx

(12)

where (a) Ni is the norm of the ith eigen function (b) ıij is the Kronecker symbol, equal to zero if i # j, or equal to 1 if i = j. By replacing Vi by its expression (8), Ni may be calculated: Ni2 =

a2i + b2i 2

+

a2i − b2i 2˛i

ai bi sin2 ˛i ˛i

sin ˛i cos ˛i +

(13)

4. Pre-treatment of the model We assume [11–13] that the solution of the previous equation system is unique but may have several equivalent mathematical expressions. We are then allowed to search one solution expression. The solution of our model is searched as a sum of two functions T+ (x,t) which includes the delayed response of the wall and T*(x,t), which is the expression of the instantaneous response of the wall. Furthermore, fast convergence of truncated expansion series is appreciated so that efficient numerical applications can be produced at the end. As the eigen functions are chosen to verify homogeneous boundary conditions, the delayed part T+ (x,t) should be constrained as well (this will avoid the Wilbraham-Gibbs phenomenon, Appendix A). Introducing the sum expression in the heat, boundary and initial Eqs. (1)–(4), we have: ∂T + ∂2 T + ∂T ∗ ∂2 T ∗ =a − + Z(x, t) +a 2 ∂t ∂t ∂x ∂x2 −

(14)

∂T ∗ Ba + Ba ∗ ∂T + T (x, t) = D(t) + T (x, t) (x, t) + (x, t) − L L ∂x ∂x x = 0 ∀t ≥ t0

+





Vi , Vj = Ni2 ıij

(7)

L

As the eigen values dimension is the inverse of a time one and as they are all negative, time constants  i of the slab are defined as the opposite-inverse of the eigen values i = −1/i = L2 /(a˛2i ). We define a scalar product in T as follows:



The orthogonality property implies:

(15)

B ∂T ∗ B ∂T + (x, t) + b T + (x, t) = E(t) − (x, t) − b T ∗ (x, t) L L ∂x ∂x x=L

∀t ≥ t0

(16)

t = t0

(17)

T + (x, t) = f (x) − T ∗ (x, t) ∀x ∈ ]0, L[

(8)

The function is chosen in order to induce homogeneous boundary conditions for T+ (x,t); the right hand side (rhs) of (15) and (16) equal to zero may be written as follows:

(9)

−

∂T ∗ Ba ∗ (x, t) + T (x, t) = D(t) L ∂x

x = 0 ∀t ≥ t0

(18)

The eigen functions are orthogonal versus the previously defined scalar product (10); they are then linearly independent and constitute a basis of T.

+

∂T ∗ B (x, t) + b T ∗ (x, t) = E(t) L ∂x

x=L

∀t ≥ t0

(19)

0

The associated norm of a function f is calculated as follows:

    f  = f1 , f2

3.1. Orthogonality of the eigen functions It can show, using (5) that if we want each eigen function verifies (1)–(3), the equation must be as follows:





Vi , Vj (i − j ) =

a [E(t)(Vj (L) − Vi (L)) + D(t)(Vj (0) − Vi (0))] L

(10)

If we want to restrict the possible expressions of the eigen functions to the orthogonal ones, the previous right hand side of Eq. (10) must be equal to zero. Among two possibilities, only one may lead to a useful way to express the eigen functions: D(t) = 0 and E(t) = 0

(11)

The eigen functions will then be chosen in order to verify individually homogeneous boundary conditions.

T* (x,t)

The boundary conditions depending on T* (x,t) are nonhomogeneous. As the initial condition (17) is used only to determine T+ (x,t), the main differential equation verified by T* (x,t), may not include any time derivative and is formed in order to simplify (14): 0=a

∂2 T ∗ + Z(x, t) ∀x ∈ ]0, L[ ∂x2

(20)

5. Instantaneous part of the solution Eq. (20) is a modified heat equation with no capacitance; this is a Poisson equation. The double integration versus x of the Poisson Eq. (20) gives the expression of the instantaneous temperature term T* : 1 T (x, t) = − a ∗

 x



Z(, t) d d + c1 (t) 0

0

x + c2 (t) ∀x ∈ ]0, L[ (21) L

2312

G. Lefebvre / Energy and Buildings 42 (2010) 2309–2322

where the integration parameters c1 and c2 are determined by introducing (21) in the boundary conditions (18) and (19):

L

c2 (t) = L

Z(, t)d +

0

Bb L

L  0

0

Z(, t)d d + a((1 + Bb )D(t) + E(t))

(22)

a (Ba + Ba Bb + Bb )

c1 (t) = −D(t)L + Ba c2 (t)

(23)

The c2 expression is based on the assumption that at least one of the Biot numbers is not equal to zero. If there is no heat source p(x, t) = 0, Z(x,t) = 0 and then: −Bb D(t) + Ba E(t) (1 + Bb )D(t) + E(t) L; c2 (t) = L c1 (t) = Ba + Ba Bb + Bb Ba + Ba Bb + Bb

(24)

and may be expressed as: Xi (t) = Xi (t0 )e

i (t−t0 )

−

ei t



t e−i 

Ni2

Vi ,

∂T ∗ ∂t

 ()d

∀t ≥ t0

(28)

t0

When the heat source term may be expressed as a separa P S (x)Uj (t)), the last two equations ble function, (T ∗ (x, t) = j=1 j become:

 dUj dXi (t) = i Xi (t) + (t) Bij dt dt P

(29)

j=1

When the heat source term Z (resp. p) may be expressed

K z (x) uk (t), (resp. p(x, t) = as a separable function Z(x, t) = k=1 k

K

p (x) uk (t)) the instantaneous term can be written versus k=1 k “natural” variables: T ∗ (x, t)



 

  

  



x x L x L x + c2Ta c1Tb + c2Tb c1ϕa + c2ϕa c1ϕb + c2ϕb ... L k

k L

  x  L  x L  1 x x 1 c1u1 + c2u1 − p1 ()d d . . . c1uK + c2uK − ... pK ()d d L L k 0 0 k 0 0 ⎡ ⎤ T1 ⎢ T2 ⎥ ⎢ ϕ1 ⎥ P ⎢ ⎥  ⎢ ⎥ Sj (x)Uj (t) × ⎢ ϕ2 ⎥ (t) = ⎢ u1 ⎥ j=1 ⎢ . ⎥ ⎣ . ⎦ . uK

=

c1Ta

(25)

with −Ba Bb Ba + Ba Bb + Bb Ba Bb = Ba + Ba Bb + Bb −Bb = Ba + Ba Bb + Bb Ba = Ba + Ba Bb + Bb

Ba (1 + Bb ) Ba + Ba Bb + Bb Bb = Ba + Ba Bb + Bb 1 + Bb = Ba + Ba Bb + Bb 1 = Ba + Ba Bb + Bb

c1Ta =

c2Ta =

c1Tb

c2Tb

c1ϕa c1ϕb



 L

L

pk ()d + Bb

L c1uk = Ba

c2ϕa

0

c2ϕb

0

k(Ba + Ba Bb + Bb )

c2uk =

 L

L



pk ()d + Bb

L

pk ()d d 0





0

pk ()d d 0

0

k(Ba + Ba Bb + Bb )

The coefficients of the temperature and flux solicitations of the instantaneous solution become simplified when the Biot numbers take their actual value (zero, non-null finite, infinite). Table 1 summarizes these expressions.

Xi (t) = Xi (t0 )ei (t−t0 ) + ei t

P 

⎛ ⎝Bij

j=1

6. Modal expansion of the delayed term The modal expansion of the delayed term is performed by searching to express it as a linear combination of the eigen functions Vi of the operator L: T + (x, t) =

∞ 

Xi (t)Vi (x)

(26)

i=1

The coefficients Xi (t) of the expansion which depend on time t are called the modal state of the thermal system. Using the orthogonal property of the eigen functions, each modal state is a solution of: 1 dXi = i Xi − 2 dt Ni



∂T ∗ , Vi ∂t



(27)

∀t ≥ t0 with Bij = −

⎞

t e−i 

dUj dt

()d ⎠

t0

1



Ni2

Sj , Vi



(30)

where U(t) is the excitation vector introduced in (25). The four first terms of the S matrix are linear versus x (Sj (x) = c1j x/L + c2j ). The projections of the unit and the x function on   an eigen function Vi are Vi , 1 = 1/˛i [ai sin ˛i − bi cos ˛i + bi ] and





Vi , x = L/˛i [(ai + bi /˛i ) sin ˛i + (ai /˛i − bi ) cos ˛i − ai /˛i ]. The four first columns of the B matrix may then calculated: Bij = −

1 Ni2 ˛2i

be

[(bi c1j + ai ˛i (c1j + c2j )) sin ˛i + (ai c1j − bi ˛i (c1j + c2j ))

cos ˛i + bi ˛i c2j − ai c1j ]

(31)

G. Lefebvre / Energy and Buildings 42 (2010) 2309–2322

2313

Table 1 Expression of sliding state coefficients versus boundary conditions (with indication of the Beck reference number). Bb = 0 Ba = 0 (X22)

No applicable

Ba = / 0 (X32)

c1 (t) = ϕb (t) +

Ba → +∞ (X12)

c1 (t) = ϕb (t) +



c2 (t) = Ta (t) +

c1uk uk (t) k 

c2 (t) = Ta (t) +

c1uk uk (t) k

/ 0 Bb = c1 (t) = −ϕa (t) +

Ba = 0 (X23)

c1 (t) =

Ba → +∞ (X13)

c1 (t) =



c2 (t) = Tb (t) +

c1uk uk (t)

Ba Bb (Tb (t)−Ta (t))−Bb ϕa (t)+Ba ϕb (t) Ba +Ba Bb +Bb Bb (Tb (t)−Ta (t))+ϕb (t) 1+Bb

c1 (t) = −ϕa (t) +

Ba = 0 (X21)

+

+





c2 (t) =

c1uk uk (t)

c1 (t) =

Ba → +∞ (X11)

c2 (t) = Ta (t) +

c1uk uk (t)

kNi2



L

Ba (Tb (t)−Ta (t))−ϕa (t) Ba +1

c1 (t) = −Ta (t) + Tb (t) +

 

ai cos

x=0

 ˛i x  L

+ bi sin

 ˛i x  L

 x 

+





c2 (t) =

c1uk uk (t)

c2 (t) = Ta (t) +

c1uk uk (t)



x

dx

0

Three integration constant series have to be determined: ai , bi and ˛i . For that purpose, we use the two boundary conditions combined with the initial condition. Introducing the general expansion expression of T+ (x,t) in the homogeneous boundary condition equations of (15) and (16), we get the following conditions which have to be verified:

ai Ba − ˛i bi = 0

1+Bb Bb

ϕa (t) +

1 Bb

ϕb (t) +

 c2uk uk (t)



+

 c2uk uk (t) k

c2uk uk (t)

 c2uk uk (t) +





c2uk uk (t) k

c2uk uk (t) k

6.1. Determination of the integration constants

sin ˛i + ˛i S cos ˛i + P sin ˛i = 0

c2uk uk (t)

Ba Ta (t)+Tb (t)+ϕa (t) Ba +1

k

(32)

−˛2i

c2uk uk (t) k

k

pj (x )dx dx 0



c2 (t) = Tb (t) + ϕa (t) +

c1uk uk (t)

Until the heat source function has to be known, only the formal expression of the other terms of the B matrix may be indicated (columns 5. . .): 1



k

k

Bij = +

ϕb (t) +

Ba (1+Bb )Ta (t)+Bb Tb (t)+(1+Bb )ϕa (t)+ϕb (t) Ba +Ba Bb +Bb

k

 k

/ 0 (X31) Ba =

1 Ba

k

k

Bb → +∞

ϕa (t) +

k

k

Ba = / 0 (X33)

1 Ba

(33) (34)

The roots are at the intersection of f(x) = B/x and g(x) = tan x curves. One root and only one exists in each range ](i − 1) , i [; the ith root will get closer to (i − 1) when i increases. When B decreases, the roots approach faster (i − 1) . The numerical root values may be calculated using for example a Newton–Raphson method or an algebraic solver such as fsolve which is available in scilab.1 ˛i cot ˛i = −B : The roots in this case are at the intersection of f(x) = −B/x and g(x) = cot x curves. One root and only one exists in each range ](i − 1) , i [; the ith root will get closer to (2i − 1) /2 when i increases. −˛2i + ˛i S cot(˛i ) + P = 0 The roots are in this case at the intersection of f(x) = (x2 − P)/xS and g(x) = cot x curves. Only one root and only one exists in each range ](i − 1) , i [; the ith root will get closer to (i − 1) when i increases. Tables 3–5 contain the first ten roots of the three previous transcendental equations for different Biot number combinations.

with S = Ba + Bb and P = Ba × Bb . The first Eq. (33) leads to many (an infinite number!) but denumerable implicit solutions ˛i : this is a transcendental equation. The second Eq. (34) allows us to determine either ai or bi in function of ˛i depending on Ba value. The remaining unknown constant is calculated using the initial condition. Depending on the Biot numbers combination, the transcendental Eqs. (33) and (34) may simplify. Table 2 resumes the solutions obtained in all but one (both Biot numbers equal to zero) configurations. Five transcendental equation types appear in Table 2. The solutions of these equations or the way to find them are now indicated.

Making the scalar product of this last equation by any eigen function and using the orthogonal property leads to:

sin ˛i = 0 : ˛i = i

Xi (t0 ) =

cos ˛i = 0 : ˛i = ˛i tan ˛i = B :

with i ∈ N

(2i − 1) 2

6.2. Initial condition The initial condition included in (4) and (17) may be written as: T (x, t0 ) =

∞ 

Xi (t0 )Vi (x) + T ∗ (x, t0 ) = f (x)

(35)

i=1

1 Ni2





(f (x) − T ∗ (x, t0 )) , Vi (x)

(36)

with i ∈ N 1 The roots values have been calculated with scilab [19]. The different codes which have been used may be downloaded from the author’s website [www.gilles.lefebvre.name].

2314

G. Lefebvre / Energy and Buildings 42 (2010) 2309–2322

Table 2 Transcendental equations to solve (“Xyz” are Beck reference numbers [8]); “?” means that the corresponding parameter value is not yet determined. Bb = 0

Bb = / 0

Ba = 0

X22

X23

Ba = / 0

X32

˛i tan ˛i = Bb

Ba → +∞

X12

ai = Bi a i and bi ? cos ˛i = 0 ai = 0 and bi ?

Bb → +∞ ˛i tan ˛i = Bb bi = 0 and ai ? −˛2i + ˛i S cot ˛i + P = 0 with S = Ba + Bb et P = Ba × Bb

X33

˛b

X21 X31

ai = ˛i bi /Ba and bi ? ˛i cot ˛i = − Bb ai = 0 and bi ?

X13

cos ˛i = 0 bi = 0 and ai ? ˛i cot ˛i = − Ba ai = ˛i bi /Ba and bi ?

X11

sin ˛i = 0 ai = 0 and bi ?

Table 3 Ten first roots of ˛i tan ˛i = B for different Biot numbers. Root index

1

2

3

4

5

6

7

8

9

10

(i − 1) B = 0.1 B = 0.5 B=1 B=2

0.000 0.311 0.653 0.860 1.077

3.142 3.173 3.292 3.426 3.644

6.283 6.299 6.362 6.437 6.578

9.425 9.435 9.477 9.529 9.630

12.57 12.57 12.61 12.64 12.72

15.71 15.71 15.74 15.77 15.83

18.85 18.85 18.88 18.90 18.95

21.99 22.00 22.01 22.04 22.08

25.13 25.14 25.15 25.17 25.21

28.27 28.28 28.29 28.31 28.34

Table 4 Roots of ˛i cot ˛i = − B for different Biot numbers. Root index

1

2

3

4

5

6

7

8

9

10

(2i − 1) /2 B = 0.1 B = 0.5 B=1 B=2

1.571 1.632 1.837 2.029 2.289

4.712 4.734 4.816 4.913 5.087

7.854 7.867 7.917 7.979 8.096

11.00 11.01 11.04 11.09 11.17

14.14 14.14 14.17 14.21 14.28

17.28 17.29 17.31 17.34 17.39

20.42 20.43 20.45 20.47 20.52

23.56 23.57 23.58 23.60 23.65

26.70 26.71 26.72 26.74 26.78

29.85 29.85 29.86 29.88 29.91

Table 5 Roots of −˛2i + ˛i S cot(˛i ) + P = 0 for different Biot numbers combinations. Root index

1

2

3

4

5

6

7

8

9

10

(i − 1) (0.5,0.1) (1.0,0.1) (0.5,0.5) (1.0,1.0)

0.000 0.734 0.929 0.960 1.307

3.142 3.321 3.453 3.431 3.673

6.283 6.377 6.452 6.438 6.585

9.425 9.488 9.540 9.530 9.632

12.57 12.61 12.65 12.65 12.72

15.71 15.75 15.78 15.77 15.83

18.85 18.88 18.91 18.90 18.96

21.99 22.02 22.04 22.04 22.08

25.13 25.16 25.18 25.17 25.21

28.27 28.30 28.31 28.31 28.35

6.3. Whole solution The following set of equations is a modal model of the heat transfer in the slab (in the particular but common case of a separable instantaneous term): T (x, t) = T + (x, t) + T ∗ (x, t)



(37)

∞

T + (x, t) =

Xi (t)Vi (x)

(38)

i=1

T ∗ (x, t) =

P 

Sj (x)Uj (t)

(39)

j=1





 Sj , Vi dUj dXi (t) = i Xi (t) − (t) dt dt Ni2 P

(2) Depending on the boundary conditions, the D(t) and E(t) functions may be evaluated. The Biot numbers may be calculated. (3) The c1 (t) and c2 (t) functions may be evaluated. (4) If the heat source function may be expressed under a separable form and the space functions are simple, the simple and double integrals versus x may be formally evaluated. If not, a numerical approximate calculation will have to be employed. (5) The instantaneous term T*(x,t) is evaluated.

(40)

j=1

As the searched solution has been split up in two terms, this model will be called a split modal model (SM). This formulation would be the same for an homogeneous slab or a whole building envelope. In the case of a slab, the evaluation will be performed more or less formally or numerically, depending on the actual form of the heat source and boundary condition functions, in the following order: 6.3.1. Instantaneous term (1) The a diffusivity is calculated.

6.3.2. Delayed term (1) The boundary conditions determine the ˛i values and possibly the ai and bi ones, or at least a relation between them. (2) The eigen function norm Ni may be calculated. (3) If the heat source function may be expressed under a separable form in which the space functions are simple, the double integral versus x may be formally evaluated. If not, a numerical approximate calculation will have to be employed. (4) The delayed term T+ (x,t) is evaluated.

6.3.3. General solution T (x, t) = T + (x, t) + T ∗ (x, t)

(41)

7. Configuration X22 (Biot numbers both equal to zero) This is the only case for which the general presented solution is not valid. As it is possible in this case to apply two incoming

G. Lefebvre / Energy and Buildings 42 (2010) 2309–2322

fluxes on the sides of the slab, the internal energy may indefinitely increase and the temperature may become (theoretically) infinite. There is no, in general, instantaneous finite term. The only particular case of no increasing internal energy is when the heat fluxes balance is equal to zero:



ϕa = −ϕb −

1

p()d

(42)

0

In general, the heat fluxes balance is positive or negative and it is impossible to formulate neither a steady-state field temperature nor a sliding state. The modal method is not applicable (as well as the variable separation method) and some special paper solutions have been put in evidence in the literature [13,14] but will not be described in this paper. 8. Dynamic simulation When the solicitations Uj (the system inputs) are known as time varying data, a dynamic simulation is the set of calculations which are done in order to know the time variations of the temperatures in the slab (outputs of the system). Among all the possible representations of the time variations, the more general and usual one consists in values taken by the variables (input solicitations or output temperatures) at regularly spaced time-steps. If t is the time-step, inputs are known and outputs are calculated at times kt for k = 1·K, K being the total number of time-steps and the initial conditions being fixed at time 0. Assuming that the (unknown) variation of the inputs Uj is linear versus time during every time-step, the integral expression of the modal states becomes: Xi ((k + 1)t)

= Xi (kt)ei t + ei (k+1)t

P 



Bij

j=1 P

= Xi (kt)ei t

Uj t



ε(x, t) =

∞ 

Xi (t)Vi (x)

(45)

i=R+1

In practice, the error value will be difficult to calculate with the previous expression as it should require knowing (a) the infinite number of the neglected eigen elements (b) and the modal state values, which depend on the way the solicitation is evoluting. Furthermore, the previous expression also depends on time t and space x. In order to eliminate the x variable, it is usual to calculate the norm of the error:

! "L ∞ ∞    " "  ε2 (t) = ε(x, t) = # Xi (t) Vi (x) Xj (t) Vj (x)dx 0 i=R+1

j=R+1

! " " ∞ 2 =# X (t)N 2 i

(46)

i

i=R+1

which has been simplified using the orthogonal property. From another point of view, we know that if we choose to observe the error during a cooling process from an initial known state with constant boundary conditions, the norm expression is:

e−i  ()d kt

(43)

1 − ei t  Uj − Bij t i

It would be possible to realize formally the integration when the sollicitations evoluate versus other simple functions such as sinusoids, polynomials (spline and Bezier functions), etc.

! " " ∞ 2 ε2 (t) = # X (0)e2i t N 2 i

(47)

i

i=R+1

As all terms of the sum are positive and the eigen values are negative numbers, the maximum error is observed at the beginning of the process:

9. Modes selection and size reduction One must pay attention to the fact that the analytical model requires extracting a finite number of terms in the infinite series of eigen elements. This modeling stage is called size reduction as it describes the dynamic thermal behavior of the slab with a reduced set of numbers. Many papers have been published about general reduction strategies of state models, which can be modal or not. Here we will mention only the simplest one, the Marshall method [15]. The eigen elements are first organized in order, i.e. the lower the index, the higher the time constant: 1 ≥ 2 ≥ · · · ≥ i−1 ≥ i ≥ i+1 ≥ · · · 1 ≥ 2 ≥ · · · ≥ i−1 ≥ i ≥ i+1 ≥ · · ·

truncature error is generated by the truncature whose expression is:

(k+1)t

j=1

0>

2315

>0

(44)

with  i = − 1/i . The Marshall method consists of a truncation of the infinite eigen elements series, keeping the R first eigen elements. When reaching steady state conditions associated with constant solicitation values U¯ j the modal states are equal to zero: Xi∞ = 0. The temperature field in steady state conditions is then expressed as:

P S (x)U¯ j . A reduction by truncation of the infinite eigen T1+ (x) = j=1 j series will have no influence on the calculation of the temperature in steady state conditions. When the thermal state is variable, a

! " " ∞ 2 ε2,max = ε2 (t = 0) = # X (0)N 2 i

i

(48)

i=R+1

Nevertheless, the initial field temperature is generally known, the initial condition being imposed to the model; we will not then observe the potential initial error, but one may be observed at long/short term. We chose to evaluate the error at 1% and 100% of the largest time constant of the slab after the initial time. When respectively t =  1 = L2 /(a 2 ) and t =  1 /100 = L2 /(100a 2 ), the errors are:

! " " ∞ 2 2 ε2 (1 ) = # Xi (0)Ni2 e−i and i=R+1

! ∞    " " 2 1 ε2 =# Xi2 (0)Ni2 e−i /100 100

(49)

i=R+1

We will see in the following example that these errors are quite small even for short expansions.

2316

G. Lefebvre / Energy and Buildings 42 (2010) 2309–2322

10. A not split modal model

and, if the stationary state is separable:

Let us apply to the previous modal equations a state variable change as follows (assuming, of course, that the solicitations may be put under a separable form):

xi (t0 ) = Xi (t0 ) −

Xi = xi +

P 

Bij Uj

(50)

j=1

P

(51)

xi (∞) = −

j=1

T + (x, t) + T ∗ (x, t)

∞ 

Xi (t)Vi (x) +

i=1 ∞

=



P 

Sj (x)Uj (t) (52)

j=1 P

xi (t)Vi (x) +

i=1

∞ 

[

j=1

Sj (x) =

Ni2

i=1

∞ 

Vi = −

Bij Vi (x) + Sj (x)]Uj (t)

10.3. Size reduction When reducing the size of a modal model by truncation, the stationary characteristics should be maintained as they determine the average behaviour and allow the model to remain energetically conservative. When reducing the modal state number of a NSM model, a static correction must be applied to the model: =

T (x, t)

i=1

T (x, t) =

Bij Vi

(53)

i=1

xi (t)Vi (x)

(54)



10.1. Initial state The initial state of the NSM verifies: 1







f (x), Vi (x) = Xi (t0 ) +

1 Ni2





T ∗ (x, t0 ), Vi (x)

(55)

R 

xi (t)Vi (x) +

∞ 

xi (t)Vi (x) +

i=1 ∞



xi (t)Vi (x)

i=R+1

(58)

xi (∞)Vi (x)

i=R+1

As it is impossible to calculate the infinite sum in the rhs of the previous equation, transformations have been made in order to get model which can be used in practice: T (x, t)



R 

xi (t)Vi (x) +

i=1 R

The model (52) and (54) is also a modal model that we will call Not Split modal model (NSM model) in order to differentiate it from the previous one. The expansion of Sj (x) in Eq. (53) is then equal to Sj (x) only over]0,L[. The idea which would consist in expanding directly the field temperature may be taken into account, but the expansion of the temperature at singular points as the slab sides may not correspond to the function (see Appendix A). As a simple illustration of the Gibbs phenomenon, in the particular X11 case, the eigen functions are sinusoids which are equal to zero on both sides of the slab. A non-null temperature may theoretically be expressed as an infinite sum of null values (the one of the eigen functions) but in practice it is impossible to reproduce this non-null temperature with a finite even large sum of null values. This explains also why a huge number of eigen functions are required to reproduce temperatures close to the sides of the slab: the convergence of the infinite series is slow. Furthermore, a huge difference will appear between both SM and NSM models (1) when expressing the initial modal state (2) when calculating the stationary thermal state (3) when applying a size reduction.

Ni2

xi (t)Vi (x) =

i=1 R

i=1

xi (t0 ) =

∞ 

i=1

and (52) becomes: ∞ 

(57)

which is generally different from zero.

Or, if Sj (x) is expansible over the {Vj (x)} eigen functions set, we have, using (29) and (40):

 ∞   Sj , Vi

Bij Uj (∞)

j=1

This equation is a simple differential equation similar to the one which determines the SM state evolutions, except that the second term of the rhs includes the sollicitations Uj but does not include its derivative. The temperature field expression becomes:

=

(56)

10.2. Stationary state

P 

 dxi (t) = i xi (t) + i Bij Uj (t) dt

=

Bij Uj (t0 )

j=1

When the slab is in a stationary state, the SM state is equal to zero and (56) indicates that the NSM state verifies:

Eq. (29) becomes:

T (x, t)

P 

=



=

xi (t)Vi (x) +

R 

i=1

 j=1

xi (t)Vi (x) +

i=1

=

xi (∞)Vi (x) −

i=1 P

i=1

R 

∞ 

P 

j=1

xi (t)Vi (x) +

P 

R 

xi (∞)Vi (x) i=1 R

Sj (x)Uj (∞) −

 Sj (x) −

 Sj (x) +

j=1

P 



i=1

R 

i=1

j=1

Sj , Vi





Ni2

Sj , Vi



Ni2

Uj (∞)Vi (x)

(59) Vi (x)

Uj (∞)

R 

Bij Vi (x)

Uj (∞)

i=1

Finally, the NSM model has been presented here only to promote the advantages of the SM model which is the one we use; among other differences which make the SM model better, it does not require any correction to be conservative and offers better convergence characteristics. 11. Validation of the SM model The validation of our solution has been performed using different methodologies; we demonstrate that the analytical expression of our solution is formally equivalent to the one of other model (analytical validation), we compare our numerical results with the ones calculated with other models (numerical validation) and we compare the numerical results produced with our model at different levels of numerical approximation (auto validation). We do not compare our numerical results with measurements (experimental validation) as it does not make sense when the validations mentioned before have been performed. Practical calculations are generally made with scilab [16].

G. Lefebvre / Energy and Buildings 42 (2010) 2309–2322

Fig. 2. Mean quadratic error due to the truncation of a modal and a Laplace solution whose infinite series have been reduced to the same order in the X11 case (center of the slab, time varying from 0 to 300 s).

Fig. 3. First eigen functions of the X11 problem.

Our method produces the following expression:

11.1. Analytical validation

11.1.1. X23 problem Özis¸ik gives the analytical solution of the X23 problem with an initial homogeneous temperature field (temperature equal to T0 ), a perfect insulation condition on the left side and a Newton condition on the right side (heat exchange coefficient hb , ambiant right temperature Tb ) obtained thanks to the separation of variables method [5, example 2-1 p. 37]. It can be shown that the analytical result given by Özis¸ik is formally equivalent to our solution: ∞   i=1

Bb

cos(˛i x/L) −a(˛ /L)2 t i e 2 2 (Bb + ˛i ) + Bb cos(˛i )



with ˛i tan ˛i = Bb



x x Ta (t) + Tb (t) L L

∞  sin(i x)

−2

i

2

Ta (0)e−a(i /L)

+2

i sin(i x) (−1) i

t 2

e−a(i /L)

+



i=1

∞ 



0 2

Tb (0)e−a(i /L)

t

dTa ()d dt

(t−)

t 2

e−a(i /L)

+ 0

i=1

=

i=0

(−1)i erfc

t

(63)

i=1

It is not obvious that both expressions are formally equivalent and we have decided to compare the numerical results that they produce. 1 calculated using an expansion limited to INF = 1000 terms is our reference ref . We calculated the quadratic mean dis1 tances E1 (R) between 1 (resp. E2 and 2 ) and ref with t varying 1 from 0 to 300 s at x = l (center of the slab) when both expansion series are limited to R elements. Fig. 2 in which results corresponding to orders 4–10 has been amplified, shows that the convergence of the modal model is better then the one of the Laplace model when reducing to orders 1,2 and 3. For orders 4, 5 and 6, the Laplace model becomes better. For orders greater than 6, quadratic mean distances are smaller than 10−3 for both models.

We calculate the SM model of a slab in the X11 case with no internal heat source: Ba → ∞ Bb → ∞

$

%

⇒

T (x = 0, t) = Ta (t) T (x = L, t) = Tb (t)

p(x, t) = 0 ∀x ∈ ]0, L[

∀t ≥ 0

(64)

∀t ≥ 0

(61) sin ˛i = 0; ai = 0; Vi = bi sin

(t−) dTb ()d dt

 i 2

√ Ni = bi 2; i = −a

L

2

i = i2 1 ; i =

˛ x i

L



;



(1 − (−1)i ) ; (i )  L 2 c with 1 = k

= i2 1 ; 1, Vi = bi

(L/˛i ) = 1 /i2 a

(65)

It is relatively easy to recognize in Fig. 3 (where bi has been taken equal to 1) every eigen function as the ith eigen function has i − 1 sign changes. The sliding term of the general solution (sliding steady state temperature field) is described by: x + c2 (t) = Sa (x)Ta (t) + Sb (x)Tb (t) ∀t ≥ 0 ∀x ∈ ]0, x[ L % x & Sa (x) = − + 1 c1 (t) = Tb (t) − Ta (t) c1a = −1 c2a = 1 L and and x c2 (t) = Ta (t) c1b = +1 c2b = 0 Sb (x) = L

T ∗ (x, t) = c1 (t)

T (x, t) − T0 1 = T∞ − T0



2

e−a(i /2l)

The only solicitations of the thermal system are temperatures Ta and Tb on respectively the left and right sides of the slab. We have:

11.1.2. X11 problem Carslaw and Jaeger [14, example 12-5, p. 308] and later Sacadura [1, p. 54] used the Laplace direct and inverse transforms in order to solve the problem of a 2l thickness wall initially at T0 and suddenly put in a fluid at T∞ . Assuming that the convective exchanges are strong, the wall is submitted to two Dirichlet conditions. If x is the space coordinate with its origin at the middle of the wall, the obtained solution expression is:

∞ 



(60)

 t

2 =

11.2. Auto validation

He gives also the analytical solution of the X11 problem with an initial homogeneous null temperature and a Dirichlet condition with time variable temperatures on both sides of the wall [5, example 5-2 p. 201] using the same method. Here again, the analytical result given by Özis¸ik is formally equivalent to our solution: T (x, t) = 1 −



 1 − (−1)i T (x, t) − T0 x =1−2 sin i T∞ − T0 i 2l ∞

Our general solution expression may be generally simplified when the actual initial and boundary conditions are known. Doing that, the simplified expression may be compared to the particular solutions proposed by other authors.

T (x, t) = 2T0

2317

&

with

 (2i + 1)l − x  √ 2 at

+ erfc

 (2i + 1)l + x  √ 2 at

(62)

(66)

2318

G. Lefebvre / Energy and Buildings 42 (2010) 2309–2322

The coefficients of the ordinary differential equations of the state space modal model are: Bia = −

2 2 =− ˛i i

Bib = +

2(−1)i i

(67)

The ordinary differential equations to be solved are then very simple: dXi 2 = i Xi − ˛i dt

 dT

a

dt

− (−1)i

dTb dt



(68)

The general temperature solution may then be expressed as: T (x, t) =

∞ 

Xi (t) sin

˛ x i

+ (Tb (t) − Ta (t))

L

i=1

x + Ta (t) L

(69)

As shown previously (43), the ordinary differential equation may be integrated in order to build a recursive relation between the modal state space value at time t + t versus the one at time t, t being the (a priori constant) time step of the dynamic simulation. The formal integration gives in the X11 case:

Xi (t + t) = Xi (t)ei t − ei t

2 ˛i



t+t

e−i 



t



dTa i dTb () − (−1) () dt dt

d

(70)

Assuming (as usual) that both solicitations Ta and Tb vary linearly versus the time t during a time step, the previous integral may be solved: Xi (t + t)

2 = Xi (Ta − (−1)i Tb ) ˛i i t = M1i Xi (t) + M2i (Ta − (−1)i Tb ) (71) (t)ei t

+ (1 − ei t )

11.2.1. Reduction by truncation Applying a Marshall reduction by keeping the R first eigen elements, the approximate field temperatures is: T˜ (x, t) =

 x

Xi (t) sin ˛i

L

i=1

x + (Tb (t) − Ta (t)) + Ta (t) L

(72)

The resulting error is: ε(x, t) = T (x, t) − T˜ (x, t) =

∞ 

Xi (t) sin ˛i

i=R+1

x L

(73)

When reaching steady state conditions associated with constant solicitations, the states and error verify: Xi (t = ∞) = 0 ⇒ ε(x) = 0

(74)

11.2.2. Case of a simple cooling process This is the X11 case in which the initial temperature field is equal to 1 and the temperature on both sides is equal to zero; we may deduct from the previous equations: T + (x) =



1−

x L



Ta +

x L

Xi (0) = 2(T0 − T + (x, t = 0)) √ 2(1 − (−1)i ) Xi (0)Ni = (i )

Tb = 0; T (x, t = 0) = 1; (1 − (−1)i ) 2(1 − (−1)i ) = ; (ibi ) (ibi )

Index

i

εR max

εR (1 /100)

εR (1 )

1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39

3.14159 9.42478 15.70796 21.99115 28.27433 34.55752 40.84070 47.12389 53.40708 59.69026 65.97345 72.25663 78.53982 84.82300 91.10619 97.38937 103.67256 109.95574 116.23893 122.52211

3.22 × 10−1 1.69 × 10−1 9.49 × 10−2 5.28 × 10−2 2.85 × 10−2 1.46 × 10−2 7.12 × 10−3 3.26 × 10−3 1.40 × 10−3 5.62 × 10−4 2.11 × 10−4 7.35 × 10−5 2.39 × 10−5 7.20 × 10−6 2.02 × 10−6 5.24 × 10−7 1.26 × 10−7 ...

3,70 × 10−5 0 ...

799 801 803 805

2510.13253 2516.41572 2522.69890 2528.98209

0.43477 0.31458 0.25795 0.22360 0.19997 0.18246 0.16880 0.15777 0.14862 0.14086 0.13418 0.12834 0.12319 0.11859 0.11445 0.11071 0.10729 0.10416 0.10128 0.09861 ... 0.01007 0.01000 0.00994 0.00988

maximum quadratic error may be calculated as:

! " INF "  2(1 − (−1)i )2

εR max = # i=R+1

(i )

2

(76)

and the quadratic error at t =  1 = L2 /(a 2 ) is:

with M1i = ei t and M2i = 2(1 − ei t )/(˛i i t).

R 

Table 6 Mean quadratic and maximum errors due to the truncation of the infinite modal serie in the X11 case (cooling from 1 to 0).

(75)

If we assume that our numerical infinite serie may be represented by a very long serie with INF terms, an approximation of the

! " INF "  (1 − (−1)i )2 e−i2

εR (t) = # 2 i=R+1

(i )

2

(77)

In Table 6, we included only the odd indices as the even ones correspond to null modal states. It can be seen that 20 states are required to get an initial quadratic error smaller than 0.1, almost 400 to reach 0.01, but only 3 states are enough to get a quadratic error at t =  1 /100 smaller than 0.01, 10 states smaller than 0.001. At time t =  1 , the quadratic error becomes negligible with only two states. Depending on the required result accuracy expressed in terms of maximum or mean quadratic error, it is possible to know the minimum reduction order to applied to the model. 11.3. An other auto-validation We simulated the dynamic behavior of a 1 m thick slab made of concrete (thermal conductivity = 2.34 W/m K, density = 2020 kg/m3 , specific heat capacity = 1000 J/kg K). The simulations have been performed over 50 time steps of 1001 s each (in case of periodic phenomena, avoiding periodicity equal to a multiple of the time step is better). The simulation duration is then equal to 13.9 h. The highest time constant is equal to 87,465 s = 24.3 h; the steady state regime will certainly not be reached during these simulations. The reference model is build with the first 1000 eigen elements when the reduced one uses only the 10 first eigen elements. Two simulation configurations are presented, both beginning with a null initial field temperature. When other values are not explicitly specified, the discussed results correspond to a model which uses 1000 eigen elements and space sampling is based on 100 nodes.

G. Lefebvre / Energy and Buildings 42 (2010) 2309–2322

2319

Fig. 4. Simulation 1 SM model (1000 eigen elements).

Fig. 6. Simulation 1; reduced SM model (10 eigen elements).

Fig. 5. Simulation 1: NSM model without static correction (1000 eigen elements).

Fig. 7. Simulation 7; reduced NSM model without static correction (10 eigen elements).

11.3.1. Simulation 1 The first simulation consists in applying a constant temperature equal to 20 ◦ C on both sides of the slab. Fig. 4 shows the field temperature (sampled over 100 nodes) at each of the 50 time steps calculated with the SM model using the first 1000 eigen elements. This will be our reference results. The same simulation performed with the NSM model (using also 1000 eigen elements) but without the static correction described by Eq. (59) illustrates the Gibbs phenomena due to the impossibility to correctly rebuild the temperature near the boundaries (Fig. 5). The evolution of the difference between the reference SM model and the NSM without static correction has been studied; the error in this case is due to the contribution of the neglected eigen elements of the field temperature; as the time constant number 1001 is equal to 87,465 s/10012 = 0.087 s, the contribution of eigen elements of order 1001 and higher, are always close to the one of the steady state regime, thus independent of time; this is the reason why only one curve has been observed and no figure is shown. Let us remind that, as the boundary temperatures are constant, these simulations are mainly a problem of initial condition memory vanishing, entirely determined by exponential terms equal to ei t = e−t/ i ; when, for example, t is greater than 5  i , the corresponding eigen element relative contribution is less than 0.7%. In our case, the first neglected eigen element has a contribution which is smaller than e−t/ i = e−1001/0.087 ! In order to point up the influence of the space sampling on the error observation, we doubled the number of nodes. Fluctuations due to high order eigen elements could be observed. The static cor-

rection is applied on the NSM model; as explained previously, the error produced by the NSM model in this case, is not numerically time dependent, and the static correction exactly compensates the error. The simulation results obtained with the NSM model are equal to those obtained with the reference model and have not been shown. 11.3.1.1. Influence of a strong reduction. A reduction at order 10 has been applied on all models. The results produced by the SM and the NSM model are so close to each other, than only the ones produced with the SM model are shown here (Fig. 6). Only the first field temperature (at time 1001 s) is a little bit different when calculated with the SM model using 1000 eigen elements or 10 elements (time constant n◦ 11 is equal to 87465 s/112 = 722 s, which means that all eigen elements of order equal or more than 11 will have a relative contribution smaller than 0.1% after 5 time steps). The NSM model which does not include any static correction (1) clearly shows that the Gibbs phenomenon is strongly disturbing the results near the sides (2) but presents huge fluctuations far from the boundaries (Fig. 7). This confirms that the NSM model must include a static correction. 11.3.2. Simulation 2 Other simulations have been performed by applying a periodically time varying temperature equal to Ta (t) = 20sin(2 t/10 × 3600) on the left side, and a constant null temperature on the right side. Fig. 8 shows the curves obtained with the reference model (SM model; 1000 eigen elements). The

2320

G. Lefebvre / Energy and Buildings 42 (2010) 2309–2322

Fig. 8. Simulation 8; SM model (1000 eigen elements).

Fig. 11. Simulation 2; reduced NSM model without static correction (10 eigen elements).

Fig. 9. Simulation 9; NSM model without static correction (1000 eigen elements). Fig. 12. Simulation 2: difference between the SM model (1000 eigen elements) and the reduced NSM model without static correction (10 eigen elements).

11.4. Numerical validation

Fig. 10. Simulation 2; difference between the SM model (1000 eigen elements) and the reduced SM model (10 eigen elements).

results obtained with the NSM model without static correction are shown in Fig. 9. The Gibbs phenomenon may be observed again near the left side of the slab. The error produced by reducing the SM model to 10 eigen elements is represented in Fig. 10 where the temperature scale has been amplified; the error is always less than 0.7 ◦ C. The same reduction applied on the NSM model without static correction produces the field temperature drawn in Fig. 11 with an error represented in Fig. 12; the largest error is 20 ◦ C on the singular point, which was known, but already 4 ◦ C at 0.1 m from the left side.

A 0.12 m thick homogeneous wall (heat conductivity = 2.39 W/m K, heat capacity = 1000 J/kg K, density = 1000 kg/m3 ) with a Newton boundary condition on its left side (heat exchange coefficient = 15 W/m2 K) and a Dirichlet one on the right side is analyzed. The temperature of the right side is maintained at 0 ◦ C when the left one varies as a pseudo random function of time (Fig. 13) taking its values in the range 0–10 ◦ C. The simulation is performed over 100 time-steps of 300 s each. The modal simulation calculates temperature fields which are sampled at 1024 equally spaced nodes; the curves produced with a 1024 modes model have been included in Fig. 14. As there is no heat source or sink, all temperatures are in the range 0–10 ◦ C. When reducing the number of kept modes, the quality of the temperature field calculation may decrease. We compare then the difference between the fields calculated with the 1024 modes model and a model whose modes number is successively 512, 256, 128, 64, 32, 16, 8, 4, 2 and 1.We calculate the mean quadratic difference at the 1024 nodes and the 100 time-steps, and the maximum absolute difference among the same space-time locations. The results are presented in Fig. 15 (where the abscissa has been limited to 16 modes as the functions do not vary significantly for higher values) where it can be seen that the error decreases quickly when the number of modes increases. The maximum quadratic error is less than 10% of the maximum amplitude and the maximum absolute error is always less then 1% (1 mode), respectively

G. Lefebvre / Energy and Buildings 42 (2010) 2309–2322

2321

Fig. 16. Mean quadratic (⊕) and maximum absolute (×) difference between results calculated with a modal model reduced at order i when using a 1024 volumes Control Volume model.

Fig. 13. Pseudo random time evolution of the left side temperature.

Fig. 14. Temperature field evolution calculated with a 1024 modes model in the X11 case with a random left side temperature evolution.

remaining difference between the results produced by both models when the number of volumes and modes increase, is due to the scheme which is used to perform the time integration. The integration of the analytical modal model is performed formally without any approximation; the sollicitations are supposed to vary linearly versus time during every time step. With the control volume model, a formal integration is not possible and an approximation is required; a purely implicite scheme is used; the results provided by the numerical model are not as good as those obtained with the analytical ones, but remains useful for some applications. Despite of the less acurate integration scheme of the numerical model, it is known to be consistant; that means that it converges toward the analytical solution when the time decreases. This can be verified when performed simulations with a 1024 volumes model sampled on 1024 nodes and a 1024 modes model using the same random excitation but with a time step divided by 5. The mean quadratic difference and the maximum absolute difference between both results are 0.0066 ◦ C and 0.0189 ◦ C, ten times smaller than in the previous case; the difference of the numerical model is then clearly explained by the time integration scheme. 12. Conclusion

Fig. 15. Mean quadratic (⊕) and maximum absolute (×) difference between results calculated with a modal model reduced at order i when using1024 modes.

A general solution of the heat conduction problem in an homogeneous plane slab has been presented; it is based on searching the solution as an expansion based on the infinite set of the eigen function of the heat operator. We showed that the more efficient model is obtained by first splitting the searched solution in a static term and a purely dynamic one. This SM model does not require any correction to be conservative and offers better convergence characteristics as it avoids the Gibbs phenomenon. We show that this analytical model may be used in practice for numerical calculations; the presented results and validation tests confirm the high quality of this model which can be considered as a reference for many others. Acknowledgements

1% and 0.1% (4 modes) and 1% and negligeable for order greater than 8 modes. This proves the fast convergence of the infinite serie. We compare now the analytical model with a numerical one based on a control volume method [17]. Volumes and nodes are chosen in order to correspond to the same points in both models. We compare analytical solutions calculated with N modes, sampled on N nodes with a control volume model using N volumes, N varying from 1024 to 1. Fig. 16 (where the abscissa is limited to 16 modes again) shows that the maximum absolute difference and the quadratic difference decrease quickly and reach an almost constant value for orders greater than 8 modes. The remaining maximum absolute difference is 0.190 ◦ C (1% of the excitation amplitude) and the quadratic one 0.034 ◦ C (0.3% of the excitation amplitude). The

The author wishes to thank Kelly Leblond and Lydie Lefebvre for their invaluable help concerning the English language. Appendix A. The Wilbraham-Gibbs phenomenon The Fourier decomposition of a function f is possible when it verifies conditions known as Dirichlet conditions but the convergence of the expansion toward the function is guaranteed only if f is a continuous function [18]. Wilbraham in 1848 [19] then Gibbs in 1899 [20] both began to comment what is now known as the Gibbs phenomenon: at discontinuity points, the Fourier expansion converges to the arithmetic

2322

G. Lefebvre / Energy and Buildings 42 (2010) 2309–2322

mean between the function values on both sides of the discontinuity. In the case of a continuous function f defined on a finite real range value [0,L] (as in our case), the Fourier expansion is made by implicitly applying an extension of the function over the whole infinite real range R in order to build a periodic function f* which corresponds to the initial function only on the definition range [0,L]. The period of extended function f* may or not be equal to the “length” L of the finite real definition range and we have no indication on how the extension is performed. It is thus not guaranteed that points x = 0 and x = L correspond to a continuity of the extended function f* . At x = 0 for example, the Fourier expansion converges toward the arithmetic mean between the unknown left limit value f(0− ) and the right limit value f(0+ ) of the extended function at the discontinuity point: lim =

f ∗ (0− ) + f ∗ (0+ ) = / f (0+ ) 2

with f (0− ) =

lim f ∗ (x) and f (0+ ) = lim f ∗ (x) = f (0). x→0 x→0 x<0 x>0

References [1] J.F. Sacadura, coord., Initiation aux transferts thermiques, Techniques and Documentation, Paris, 1980 (Chapter 2). [2] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, second ed., Clarendon Press, Oxford, 1959 (Chapter III). [3] M.N. Özis¸ik, Boundary Value Problems of Heat Conduction, Dover Pub, New York, 1989.

[4] M.D. Mikhailov, M.N. Özis¸ik, Unified Analysis and Solutions of Heat and Mass Diffusion, Dover Pub, New York, 1993 (§4.1 and Chapter 7). [5] M.N. Özis¸ik, Heat Conduction, John Wiley & Sons, New York, 1980 (Chapters 2, 5–7). [6] J. Martinet, Thermocinétique approfondie, Techniques et Documentation, Lavoisier ed, Paris, 1990, Chapter 1–5. [7] F. de Monte, Multi-layer transient heat conduction using transition time scales, International Journal of Thermal Sciences 45 (September (9)) (2006) 882–892. [8] G. Pontrelli, F. de Monte, Mass diffusion through two-layer porous media: an application to the drug-eluting stent, International Journal of Heat Mass Transfer 50 (August (17–18)) (2007) 3658–3669. [9] J.V. Beck, et al., Heat Conduction Using Green’s Functions, Series in Computational Methods in Mechanics and Thermal Sciences, Hemisphere Pub, 1992 (Chapter 2). [10] K.D. Cole, Organize Green Functions, http://www.engr.unl.edu/∼glibrary, September 2007. [11] M. Fourier, Théorie (analytique) de la chaleur, Chez Firmin Didot Père et Fils, Paris, 1822. [12] G. Lefebvre, La méthode modale en thermique, Technosup, Ellipses éd, Paris, 2007, idem (Annex A). [13] J. Génot, Thermique in Encyclopédie Universalis, 9, Paris, 1985, pp.1145–1154. [14] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, Oxford at the Clarendon Press, 1947. [15] S.A. Marshall, An approximate method for reducing the order of a linear system, Control 10 (1966) 642–643. [16] http://www.scilab.org, 2010. [17] G. Lefebvre, Modal-based simulation of the thermal behavior of a building: the m2m software, Energy and Buildings 25 (January (1)) (1997) 19–30. [18] R. Radaelli-Sanchez, Dirichlet conditions, v2 9, http://www.cnx.org/content/ m10089/latest/may, 2008. [19] H. Wilbraham, On a certain periodic function, Cambridge & Dublin Math. J., 3, 1848, p. 198. [20] J.W. Gibbs, Fourier’s series, Nature 1539 (59) (1899) 606.