Derivation of the heat transfer matrix inside ISO 13786

Spread the love

The purpose of this article is to provide the derivation of the heat transfer matrix of ISO 13786 from the solution of the one-dimensional heat equation (partial differential equation, PDE). The matrix refers to building components as defined in ISO 13786. Inside the same regulation the heat transfer matrix is given without a proof and no derivation is available online.

In the one-dimensional heat equation considered in this case, T(x,t) is the temperature which is a function of the x-coordinate and time.

(1)    \begin{equation*}\frac{\partial T}{\partial t} = \frac{\lambda}{c\rho} \frac{ \partial^2 T}{\partial x^2} \end{equation*}

This partial differential equation is solved by the so-called method of separation of variables. The solution is supposed to be made of the product of two other functions of one of the variables x and t separately:

(2)    \begin{equation*} T(x,t) = T_0 + U(t) \cdot X(x) \end{equation*}

The constant inside the T(x, t) function is neglected:

(3)    \begin{equation*} \Delta T(x,t) = T(x,t) - T_0 = U(t) \cdot X(x) \end{equation*}

The PDE becomes then:

(4)    \begin{equation*}\frac{\partial \Delta T}{\partial t} = \frac{\lambda}{c\rho} \frac{ \partial^2 \Delta T}{\partial x^2} \end{equation*}

Now, the T(x,t) function is replaced by the two separate functions inside the PDE and the product rule is employed:

(5)    \begin{equation*} X(x) \cdot U'(t) =  \frac{\lambda}{c\rho}  U(t) \cdot X''(x) \end{equation*}

The above equation is rearranged as follows:

(6)    \begin{equation*} \frac{X''(x)}{X(x)} =  \frac{c\rho}{\lambda}  \frac{U'(t)}{U(t)} = \mu \end{equation*}

We assume now that μ is a complex number with the real part equal to 0.

(7)    \begin{equation*} \mu \in, \mathbb{C},  \Re(\mu) = 0 \end{equation*}

Now the above PDE can be split into two separate ordinary differential equations (ODE) which can be solved separately:

(8)    \begin{equation*}  \begin{cases} X''(x) = \mu X(x) \\ U'(t) = \frac{\lambda \mu}{c \rho}  U(t)\end{cases}\end{equation*}

(9)    \begin{equation*}  \begin{cases} X(x) = A e^{\sqrt{\mu} x} + A e^{-\sqrt{\mu} x} \\ U(t) = B e^{\frac{\lambda \mu}{c \rho} t}\end{cases}\end{equation*}

(10)    \begin{equation*} \mu = k \cdot j \end{equation*}

(11)    \begin{equation*}  \begin{cases} X(x) = A e^{\sqrt{k j} x} + B e^{-\sqrt{k j} x} \\ U(t) = C e^{\frac{\lambda k j}{c \rho} t} = B e^{j \omega t}\end{cases}\end{equation*}

We introduce the new variable ω (angular velocity), which is equal to:

(12)    \begin{equation*} \omega =  \frac{\lambda k}{c \rho} \implies k = \frac{\omega c \rho}{\lambda} \end{equation*}

The power of the exponentials inside the X(x) function can be rewritten as:

(13)    \begin{equation*} \sqrt{k j} = \sqrt{\frac{\omega c \rho}{\lambda} j} = \sqrt{\frac{\omega c \rho}{\lambda} e^{j \frac{\pi}{2}}} = \sqrt{\frac{\omega c \rho}{\lambda}} e^{j \frac{\pi}{4}} =  \sqrt{\frac{\omega c \rho}{\lambda}} \left ( cos \left (\frac{\pi}{4} \right ) + j sin \left ( \frac{\pi}{4} \right ) \right ) = \end{equation*}

(14)    \begin{equation*} = \sqrt{\frac{\omega c \rho}{\lambda}} \left ( \frac{1}{\sqrt{2}} + j \frac{1}{\sqrt{2}} \right ) = \sqrt{\frac{\omega c \rho}{2\lambda}} ( 1 + j) = \frac{1}{\delta} (1+j) \end{equation*}

δ is defined as:

(15)   \begin{equation*}\delta =  \sqrt{\frac{2\lambda}{\omega c \rho}}\end{equation*}

The two differential equations above can be written as:

(16)    \begin{equation*}  \begin{cases} X(x) = A e^{\frac{1}{\delta} (1+j) x} + B e^{-\frac{1}{\delta} (1+j) x} \\ U(t) = C e^{j \omega t}\end{cases}\end{equation*}

Now, the solution of the one-dimensional heat equation can be written:

(17)    \begin{equation*}  \Delta T(x,t) =\left (  A e^{\frac{1}{\delta} (1+j) x} + B e^{-\frac{1}{\delta} (1+j) x} \right )  C e^{j \omega t} = \left (  A' e^{\frac{1}{\delta} (1+j) x} + B' e^{-\frac{1}{\delta} (1+j) x} \right )  e^{j \omega t} \end{equation*}

where

(18)    \begin{equation*}  A' = A \cdot C\end{equation*}\begin{equation*}  B' = B \cdot C\end{equation*}

The following parameters are introduced:

θ0: maximum temperature at side n of the building component
q0: maximum heat flux density at side n of the building component
θL: maximum temperature at side m of the building component
qL: maximum heat flux density at side m of the building component

The origin of the x-axis is located on side n while the x-value on the opposite side m is equal to L. Therefore, the solution can be written on the two boundaries as:

(19)    \begin{equation*} \Delta T(0,t) = \theta_0 e^{j \omega t} = (A' + B') e^{j \omega t} \implies \end{equation*}

(20)    \begin{equation*} \implies \theta_0 = (A' + B') \end{equation*}

(21)    \begin{equation*} q(0,t) = q_0 e^{j \omega t} = -\lambda \frac{\partial T}{\partial x}  \bigg|_{x=0} = -\lambda \left { \left [ A' \frac{1}{\delta} (1+j) e^{\frac{1}{\delta} (1+j) x} - B' \frac{1}{\delta} (1+j) e^{- \frac{1}{\delta} (1+j) x} \right ] e^{j \omega t} \right } \bigg|_{x=0} = \end{equation*}

(22)    \begin{equation*} = - \frac{\lambda}{\delta} (1+j) (A'-B') e^{j \omega t} \implies q_0 = \frac{\lambda}{\delta} (1+j) (A'-B') \end{equation*}

(23)    \begin{equation*}  \Delta T(L,t) = \theta_L e^{j \omega t} = \left [  A' e^{\frac{1}{\delta} (1+j) L} + B' e^{-\frac{1}{\delta} (1+j) L} \right ]  e^{j \omega t} \implies \end{equation*}

(24)    \begin{equation*} \implies \theta_L = \left [  A' e^{\frac{1}{\delta} (1+j) L} + B' e^{-\frac{1}{\delta} (1+j) L} \right ] \end{equation*}

(25)    \begin{equation*} q(L,t) = q_{L} e^{j \omega t} = -\lambda \frac{\partial T}{\partial x}  \bigg|_{x=L} = -\lambda \left { \left [ A' \frac{1}{\delta} (1+j) e^{\frac{1}{\delta} (1+j) x} - B' \frac{1}{\delta} (1+j) e^{- \frac{1}{\delta} (1+j) x} \right ] e^{j \omega t} \right } \bigg|_{x=L} \end{equation*}

(26)    \begin{equation*} \implies q_{L} = -\lambda \left [ A' \frac{1}{\delta} (1+j) e^{\frac{1}{\delta} (1+j) L} - B' \frac{1}{\delta} (1+j) e^{- \frac{1}{\delta} (1+j) L} \right ] \end{equation*}

The four equations obtained are as follows:

(27)    \begin{equation*}  \begin{cases} \theta_0 = (A' + B') \\ q_0 = \frac{\lambda}{\delta} (1+j) (A'-B') \\  \theta_L = \left [  A' e^{\frac{1}{\delta} (1+j) L} + B' e^{-\frac{1}{\delta} (1+j) L} \right ] \\ q_{L} = -\lambda \left [ A' \frac{1}{\delta} (1+j) e^{\frac{1}{\delta} (1+j) L} - B' \frac{1}{\delta} (1+j) e^{- \frac{1}{\delta} (1+j) L} \right ] \end{cases}\end{equation*}

Now we rearrange the first two equations so that parameters A’ and B’ are isolated on left-hand side of each of the two equations:

(28)    \begin{equation*}  \begin{cases} A' = \frac{\theta_0}{2} - \frac{q_0}{4\lambda} \delta (1-j) \\  B' = \frac{\theta_0}{2} + \frac{q_0}{4\lambda} \delta (1-j) \end{cases}\end{equation*}

Those two parameters are then replaced in the third and fourth equations:

(29)    \begin{equation*}  \begin{cases} \theta_L = cosh \left ( \frac{L}{\delta} (1+j)  \right ) \theta_0 - \frac{\delta}{2 \lambda} (1-j) sinh \left ( \frac{L}{\delta} (1+j)  \right ) q_0 \\ q_L = - \frac{\lambda}{\delta} (1+j) sinh \left ( \frac{L}{\delta} (1+j) \right ) \theta_0 + cosh \left ( \frac{L}{\delta} (1+j) \right ) q_0 \end{cases} \end{equation*}

The following properties of hyperbolic functions with complex argument are used in order to obtain the final form of the matrix provided inside regulation ISO 13786:

(30)    \begin{equation*}  sinh(z_1 + z_2) = sinh(z_1)cosh(z_2) + cosh(z_1)sinh(z_2) \end{equation*} \begin{equation*}  cosh(z_1 + z_2) = cosh(z_1)cosh(z_2) + sinh(z_1)sinh(z_2) \end{equation*} \begin{equation*}  cosh(ix) = cos(x) \end{equation*} \begin{equation*}  sinh(ix) = i \cdot sin(x) \end{equation*}

After a few calculations and after relacing parameter ξ, defined as:

(31)    \begin{equation*}  \xi = \frac{L}{\delta} \end{equation*}

the final form is obtained:

(32)   \begin{equation*} \left [ \begin{array}{c} \theta_L \\ q_L \end{array} =  \left [ \begin{array}{cc} M_{1,1} & M_{1,2} \\ M_{2,1} & M_{2,2} \end{array} \right ] \cdot \left [ \begin{array}{c} \theta_0 \\ q_0 \end{array} \right ] \end{equation*}

where:

(33)    \begin{equation*}  M_{1,1} = M_{2,2} = cosh(\xi)cos(\xi) + j sinh(\xi) sin(\xi) \end{equation*} \begin{equation*}  M_{1,2} = \frac{\delta}{2\lambda} \left [ sinh(\xi)cos(\xi) + cosh(\xi)sin(\xi) + j \left ( cosh(\xi) sin (\xi) - sinh(\xi) cos(\xi) \right ) \right ] \end{equation*} \begin{equation*}  M_{2,1} = - \frac{\lambda}{\delta} \left [ sinh(\xi)cos(\xi) - cosh(\xi)sin(\xi) + j \left ( cosh(\xi) sin (\xi) + sinh(\xi) cos(\xi) \right ) \right ] \end{equation*}

Relationship between complex solution and spatial solution

References

Parametri Dinamici UNI EN ISO 13786 – Prof. Marco Manzan

Luca Bonifazi – Comportamento termico dinamico di strutture opache