Thermal conduction in dynamic conditions

The problem of conduction heat flow can be solved by means of direct integration of the Fourier differential equation. As an example, let us consider a semi-infinite plane, with a boundary condition of a sinusoidal temperature of t = t₀ cos ωτ at one face (Figure 3.8.a).

If the axe x has its origin in this surface and direction in the verse to the body, the equations of temperature and heat flux (3.11) and (3.12) can be written as:

(

x

)

e t

t = ₀ ⋅ ^−βx cosωτ −β (3.26)

(

^/4

)

cos

2λβ ₀ ^β ωτ −β +π

= t e⁻ x

q ^x (3.27)

where

a 2

β = ω . From equation (3.26) the temperature inside the body is not in phase, with a delay time of βx/ω (depending on the deepness) and the amplitude decreases in the deepness with an exponential law (e^-βx). Equation (3.27) shows that the flux has an anticipation of a quarter of period on the temperature and an amplitude proportional to the temperature by means of a multiplication factor 2λβ.

The presence of discontinuity in the material produces reflection waves which complicate the analytical solution of equations (3.11) and (3.12).

56

a b

Figure 3.8: Problem of sinusoidal boundary condition(a) and temperature trend in the thickness of the semi-infinite medium (b)

A solution of this problem is to use the electrical analogy. The thermal problem can be modelled by means of an electrical circuit [27, 28] with series of resistances and capacities. As a consequence of the slowness of the thermal phenomenon (which has a delay as a function of the co-ordinate even with low values of thickness) the electrical circuit has to be defined by means of distributed constants, i.e. the circuit is characterised by a specific resistance r and capacity c. In the conductor the self-induction and the dispersion conductance can be considered negligible. The equations for a such system can be written as:

x

Let us consider again a surface temperature (electromotive force in the electrical circuit) with a sinusoidal law. This case can be solved by means of Fourier series, resorting to the symbolic method; in this way equations can be written in the vector way:

x q

These equations can be written as:

t

57

attenuation constant and β the phase constant which defines the periodicity in the space domain; Z& =1 jωλρc_p is the impedance characteristic of a thermal conductor.

A multi-layer wall can be modelled by means of a series of electrical circuits (Figure 3.9).

In this way, if V&_i,V&_o,I&_i,I&_o are respectively incoming and outgoing voltage and the intensity current, it follows that:



with A&,B&,C&,D& constants which depend on the elements constituting the electrical circuit and on their distribution. Considering a simple symmetric circuit with T shape, due to Kirchhoff principles the constants can be expressed as:

2 ,

This circuit can be considered equivalent to a heat conductor with finite length l in which the impedances can be written as:

l R

Figure 3.9: Electrical analogy for a multi-layer wall

Considering more circuits in series, each of them characterised by the thermal impedances characteristic of the material, temperature and heat flux in one surface can be evaluated by means of equation (3.36), if the other boundary condition changes with a sinusoidal law. The wall can be solved if the constants of the whole system is calculated as the product of the matrices of the constants of each layer. The product has to start from the layer whose boundary conditions are known; in this way a unique equivalent circuit in series is possible to define. Similar considerations can be pointed out for a circuit in parallel, but this case is not very interesting for the building applications.

For each system, knowing two of the functions t&_i,q&_i,t&_o,q&_o the other two can be evaluated by means of the vector representation, as attenuation and delay in phase with respect to the boundary conditions imposed.

58 3.3.2 Numerical solutions

Numerical methods can be a useful tool to solve the problems of thermal conduction in transient conditions. The most important are the finite difference, the finite element and the heat transfer function methods.

3.3.2.1 Finite difference method

The finite difference method is based on the substitution of the differential equation (3.9) with equations with finite differences. Let us supposed to divide the whole thickness of the medium in layers of thickness ∆x (each space is called “node”). At a given time τ the temperatures of three generic adjacent nodes is t_i-1, t_i, t_i+1. At time step τ+∆τ the temperatures are t’i-1, t’i, t’i+1. Considering the heat exchange located in these adjacent nodes and making the assumption that the thermal characteristics of the layers (with thickness ∆x) are concentrated in the nodes, the thermal balance of the i-th node in the time range between τ and τ+∆τ can be written in these two ways:

Introducing the for time and space steps the above equations can be written as:

(

_{i i i}

)

which can be used respectively if the temperatures in the second member at time τ are calculated in the previous step or if the ones at time τ+∆τ have to be calculated in the same step (and therefore are still unknown). The equation (3.38) is called explicit method, the equation (3.39) is called implicit method.

Figure 3.10: Scheme of finite difference technique in one dimension problems For each of the n nodes it is necessary to write the equation (3.38) or (3.39), depending on the method adopted; in this way a system with n equations and n variables is obtained. In the first case, there is one variable for each equation, so the equations can be solved alone. In the second method there is the need to solve the problem of the whole system, but it has the main advantage to have a matrix three-diagonal, which means that some calculation procedures “ad hoc” can be assumed.

59

Since the original differential equation is substituted by approximating equations, there are some truncation errors. A possible way to evaluate them is by means of the Taylor series to be compared with the exact equations. It is possible to demonstrate that these errors are in the order of ∆τ for the time variable and in the order of ∆x² for the space variable (i.e. order 1 in ∆τ and order 2 in ∆x).

As far as the numerical stability of the methods is concerned, the explicit method has a limitation in the choice of the variables. In fact, making the assumption that at time step τ the heat flux in the i-th node is positive, i.e. the left side of equation (3.38), it means that: again positive or zero:

0 stability criterion for choosing the correct ∆x e ∆τ . In the implicit method temperatures are calculated at the end of the time step; therefore flux inversions are not possible and therefore the second principle is always satisfied.

There are some other one step methods, implicit or explicit, like the ones already shown.

For example Crank-Nicholson method is implicit of the second order with respect to both the variables.

3.3.2.2 Finite element method

This method is based on the subdivision of the system in a finite group of element, by recurring to an interpolation analytical function. This method is different from the finite difference technique as the use of matrices is always necessary.

Let us consider equation (3.11). The boundary conditions for the problem are:

t1

where Sa and Sb are the boundary surfaces of the one-dimensional system.

Galerkin’s approach is one of the most common. The solution domain is divided into finite elements in space (Figure 3.11). The temperature is approximated within each element nodal temperature considered to be dependent on time and m is the number of nodes in the considered element.

60

Figure 3.11: Scheme of finite element technique in one dimension problems The Galerik representation for the heat conduction problem in equation (3.11) is

=0

Using integration by parts in the first three terms, equation (3.44) simplifies to

=0

Inserting the temperature approximation, equation (3.45) becomes

{ } { }

The above equation can be cast into a more convenient form as

{ }

[ ]

{ }

0

This set of equation requires an iterative solution. Following the simplest form of iteration method some initial attempt has to be assumed:

(

₁⁰ ₂^{0 0}

)

0 t ,t ,...,t_m t

t = = (3.50)

and obtain an improved solution t¹ by solving the equation

{ }

^{[ ]}

{ }{ }

⁰

The general iteration scheme

61

{ }

^{[ ]}

{ }{ }

⁰

]

[ ⁿ + K tⁿ⁻¹ tⁿ = d

t M d

τ ^(3.52)

is then repeated until convergence, within a suitable tolerance, is obtained.

3.3.2.3 Transfer functions method

Response factors are based on transfer functions D, which relates the input Ω (τ) applied to a certain physical system to the response O(τ). The law relating the input and the response is [30]:

) (

* )

(τ = D Ωτ

O (3.53)

where * is a convolution in the time domain.

The functions O(τ) and Ω (τ) are usually continuous in time and can be substituted by the related series. The series of a function F(τ), which is valid between τ 0 and τ 0 + b ∆τ , is a set of the b terms Fj (j = 1, … b) obtained by sampling the values of the function in a time range ∆τ. Considering the triangle with height equal to the value of the function and base equal to two times the time step (Figure 3.12), the superposition of these triangles can give a satisfactory approximation of the original function F(τ).

Equation (3.53) allows to define the transfer function D as the response of the physical system due to an impulse Ω u(τ), which corresponds to a triangle with a unitary height. In general D can be represented by a temporal series D_j (j = 1, … ∞ ) with an infinite number of terms, as the response of a system to an impulse is asymptotic (Figure 3.13).

Figure 3.12: Approximatio of the function F(τ) by means

of the series of triangles Fj

Source: [30]

Figure 3.13: Example of shape of the response function Dj due to the triangular impulse Ω u

Source: [30]

As the generic term Dq of the transfer function represents the response of the system to a unitary triangular impulse applied in an instant before q ∆τ the considered one.

Superposing the effects of each term of the effective impulse, the value of the response Ok at the instant τ = k ∆τ can be expressed as:

∑

∞

= Ω − +

=

, 1

1 j

j k j

k D

O (3.54)

The elements Dj (j = 1, …

∞

) of the transfer function are called “response factors”. In practice it is recommended that the series (3.54) would be stopped at a finite number of terms; in order to achieve a good precision of the method, a sequence of modified terms

62

should be hence evaluated. In the succession of response factors (whose trend is asymptotic to zero), it is always possible to evaluate a number N, so that for each term greater than N the ratio between two consecutive terms can be considered constant, i.e.

j R

j c

D

D ⁺¹ ≈ , where cR is a constant called “common ratio”.

It can be demonstrated that equation (3.54) can be substituted with:

=

∑

Ω − + + −

=

N j

k R j

k j

k D c O

O

, 1

1

' 1 (3.55)

where D'_j are the modified factors of the transfer function, whose values are:

1

'1 D

D = (3.56a)

'_j =D_j −c_RD_j−1

D (3.56b)

Equation (3.55) allows to evaluate Ok by means of a finite number of known terms.

In problems related to heat transfer the use of transfer function deals with heat fluxes and temperatures, where the temperature is the input and the heat flux is the response of the system.

In the most simple case (1-D problem), only two boundary conditions are present: the internal surface temperature and the outer temperature (usually, the external sol-air or the adjacent ambient temperature). In this case, the heat flux is actually one dimensional and the transfer functions required will be defined (and obtained) as follows (see Figure 3.14):

Function Z_i→i: response of conductive heat flux at the inner surface of the structure, due to an unitary triangular impulse of temperature on the same surface, keeping the outer temperature constantly equal to zero.

Function Z_o→i: response of conductive heat flux at the inner surface of the structure, due to an unitary triangular impulse of temperature on the outer surface, keeping the internal temperature constantly equal to zero.

Function Zi→o: response of conductive heat flux at the outer surface of the structure, due to an unitary triangular impulse of temperature on the inner surface, keeping the outer temperature constantly equal to zero.

Function Zo→o: response of conductive heat flux at the outer surface of the structure, due to an unitary triangular impulse of temperature on the outer surface, keeping the inner temperature constantly equal to zero.

It is easy to show that, for symmetry reasons, the transfer function Zo→i is equal to Zi→o.

Figure 3.14: Transfer functions for building elements not containing pipes Impulse

Response Z Response Z

o i

o o

- - ^{Response Z}

^{i o}

^{Response Z}

^{i i}

Im pulse

-

-63

The equation of the heat flow for the inner surface q*i and for the outer surface q*o, with surface temperatures t_si, and t_so respectively, can be written in the time step τ=p∆τ: leading to the equations:

( ) ( ) ( ) ( )

₁

The transfer functions method can be applied by means of the Z-transform functions to solve one dimension conduction thermal flow. In this way, the calculated response factors can be implemented for heat transfer across walls in the thermal balance of a room. This calculation procedure has been developed by Mitalas, Kusuda and Stephenson [31, 32, 33, 34, 35] via Laplace transform functions. This technique sometimes can lead to no satisfactory results. For example when the structure is well insulated or thick, or if the time step for simulations is very short. The first aspect is very important if a study has to be lead on the effect of increasing insulation in buildings or in historical buildings where the walls can be 1 m or 2 m thick. The second aspect is important as well, because when looking at the control strategy of a HVAC for heating/cooling a building, time steps shorter than 1h are needed to simulate more accurately the transient conditions of the plant in the rooms.

As an alternative, thermal response factors can be also evaluated by means of the finite difference or the finite elements technique, which require longest calculation time, but are accurate and allow to solve the problem even with short time steps or with thick walls.

3.2.3 Building elements containing pipes

As already explained, when considering a radiant system with pipes embedded in the structure it is common to assume that the middle section between the two pipes is adiabatic (Figure 3.7). This hypothesis is acceptable, since the distance between the two pipes is generally between 10 and 20 cm. Thus, the heat flow between the pipes is negligible with respect to the heat flow through the section.

One possible way to simulate heat transfer in dynamic conditions within the structure of a radiant system is to use the response factor technique. In this case three simulations are necessary to fully describe the thermal conduction in a building structure with embedded pipes: in each simulation the triangular impulse is given on the inner surface, the outer surface, or on the internal surface of the pipes. As a result of the simulation, the heat flow on each surface is recorded (Figure 3.15).

64

Figure 3.15: The transfer functions inside a structure that contains embedded pipes The behaviour of a building element with embedded pipes can be therefore expressed by means of nine response factors:

• Function Z_i→i: the thermal response of the specific flux on the internal surface (i) due to a unitary impulse in the temperature on the same surface (i), which keeps the temperature on the external surface and on the pipe’s internal surface constantly zero.

• Function Z_o→i: the thermal response of the specific flux on the internal surface (i) due to a unitary impulse in the temperature on the outer surface (o), which keeps the temperature on the internal surface and on the pipe’s internal surface constantly zero.

• Function Z_p→i: the thermal response of the specific flux on the internal surface (i) due to a unitary impulse in the temperature on the pipe’s internal surface (p), which keeps the temperature on the internal and on the external surfaces constantly zero.

• Function Z_i→o: the thermal response of the specific flux on the outer surface (o) due keeps the temperature on the internal and on the external surfaces constantly zero.

• Function Z_i→p: the thermal response of the specific flux on the pipe’s internal which keeps the temperature on the internal and on the external surfaces constantly zero.

65

calculated by means of a commercial software programme [36] that is based on the FDM.

It may be seen that Z_i→o = Z_o→i, Z_i→p = Z_p→i, Z_o→p = Z_p→o, therefore the functions that need to be calculated are six instead of nine.

By superimposing the effects, it is possible to obtain the specific heat fluxes on the internal surface (q*_i), the outer surface (q*_o) and the internal surface of the pipe (q*_p), for a generic series of temperatures on the internal surface (ts,i), on the outer surface (ts,o) and on the internal surface of the pipe (ts,p).

Once the six response factors have been obtained, the common ratio, cR, must be calculated in order to evaluate the modified response factors:

( )

factors for the impulse on the pipe surface are reported; it can be observed that in this case 71 modified response factors are necessary to describe the tile covered radiant floor [37].

In this way, the heat fluxes due to conduction through a structure with embedded pipes (q*i, q*o and q*p) can be written via the modified response factors in the following way:

66 a)

b)

c)

Figure 3.16: Response factors for tile covered radiant floor: from internal to outer surface (a), from internal to pipe surface (b), from pipe to outer surface (c)

-0.020 -0.018 -0.016 -0.014 -0.012 -0.010 -0.008 -0.006 -0.004 -0.002 0.000

0 10 20 30 40 50 60 70 80 90 100

[W/(m2 K)]

Z i-o Z o-i

-7 -6 -5 -4 -3 -2 -1 0

0 10 20 30 40 50 60 70 80 90 100

[W/(m2 K)]

Z i-p Z p-i

-0.020 -0.018 -0.016 -0.014 -0.012 -0.010 -0.008 -0.006 -0.004 -0.002 0.000

0 10 20 30 40 50 60 70 80 90 100

[W/(m2 K)]

Z p-o Z o-p

67 a)

b)

c)

Figure 3.17: Example of response factors and modified response factors obtained for a typical tile covered radiant floor: from pipe to pipe surface (a), from internal to pipe

surface (b), from pipe to outer surface (c)

-20 -15 -10 -5 0 5 10 15 20

0 10 20 30 40 50 60 70 80

[W/(m2 K)]

Z p-p Z' p-p

-7 -6 -5 -4 -3 -2 -1 0 1 2

0 10 20 30 40 50 60 70 80

[W/(m2 K)]

Z i-p = Z p-i Z' i-p = Z' p-i

-0.020 -0.018 -0.016 -0.014 -0.012 -0.010 -0.008 -0.006 -0.004 -0.002 0.000

0 10 20 30 40 50 60 70 80

[W/(m2 K)]

Z p-o = Z o-p Z' o-p = Z' p-o

68

An example to show the accuracy of the response factor method is reported in Figure 3.18, where the results of a simulation are shown in terms of heat flow through the pipe surface [37]. In this simulation the pipe surface temperature changes from 20°C to 70°C (step-wise constant), while the internal and external surfaces remain constantly at 20°C.

As shown, the results of the response factor method are the same of the detailed FDM, therefore the response factor method can accurately predict the heat conduction through building elements containing pipes.

Figure 3.18: Example of comparison between FDM and response factor technique in terms of heat flow through the pipe surface [W/m²] due to 50°C sudden change in the

surface temperature of the pipe (the inner and outer surface remain at constant temperature)

Comparison between detailed calculations (FDM) and response factors technique

-400 -200 0 200 400 600 800 1000

0 250 500 750 1000 1250 1500 1750 2000

Minutes Heat flow [W/m2 ]

20 30 40 50 60 70 80 90

Temperature [°C]

FDM Response Factors Temperature in the pipe

69

In document Simulation and numerical methods Energy modeling for buildings and components (Pldal 55-69)