Home Heating

Home → Home Heating

Any owner will confirm that the heating of a house is expensive. In addition to the cost, it's also a carbon pollution! It turns out that using a system of differential equations, it is possible to describe the change of the temperature in the house for different heating modes. Such a model can be of real practical value.

Consider a typical two-storey house (Figure $1$ ).

We will describe the indoor temperature by two functions:

$x (t)$ is the temperature of the ground floor;
$y (t)$ is the temperature of the 1st (upper) floor.

The initial differential equations are based on the Newton's law of cooling, which is written as

\frac{d T}{d t} = \frac{α A}{C} (T_{e} - T) = k (T_{e} - T),

where $T_{e}$ is the temperature of the surrounding environment or adjacent room, $T$ is the temperature of the room, $α$ is the heat transfer coefficient, $C$ is the heat capacity of the body or surface through which the heat passes, $A$ is the surface area. It is convenient to replace the coefficients $α,$ $A$ and $C$ with one parameter − the coefficient of thermal conductivity $k,$ which describes the rate of temperature equalization between the two points. Thus, when $k = 1,$ the initial temperature difference $Δ T$ decreases $e \approx 2.7$ times per unit time. With increasing $k,$ the relaxation rate also increases.

We introduce the following coefficients $k_{i}$ describing the two-storey house:

$k_{1}$ is the thermal conductivity of the floor on the ground floor;
$k_{2}$ is the thermal conductivity of the ceiling on the ground floor;
$k_{3}$ is the thermal conductivity of the walls on the ground floor;
$k_{4}$ is the thermal conductivity of the walls on the first floor and roof.

Then the system of differential equations for the functions $x (t),$ $y (t)$ can be represented as follows:

{\begin{cases} \frac{d x}{d t} = k_{1} (T_{g} - x) + k_{2} (y - x) + k_{3} (T_{e} (t) - x) + f (t) \\ \frac{d y}{d t} = k_{4} (T_{e} (t) - y) + k_{2} (x - y) \end{cases}, \Rightarrow {\begin{cases} \frac{d x}{d t} = - (k_{1} + k_{2} + k_{3}) x + k_{2} y + k_{1} T_{g} + k_{3} T_{e} (t) + f (t) \\ \frac{d y}{d t} = k_{2} x - (k_{2} + k_{4}) y + k_{4} T_{e} (t) \end{cases},

where the term $f (t)$ describes the heat source located on the ground floor.

We write this system in vector-matrix form:

X^{'} (t) = K X (t) + F (t), where X (t) = [\begin{matrix} x (t) \\ y (t) \end{matrix}], K = [\begin{matrix} - k_{1} - k_{2} - k_{3} & k_{2} \\ k_{2} & - k_{2} - k_{4} \end{matrix}], F (t) = F_{1} (t) + F_{2} (t) = [\begin{matrix} k_{1} T_{g} + k_{3} T_{e} (t) \\ k_{4} T_{e} (t) \end{matrix}] + [\begin{matrix} f (t) \\ 0 \end{matrix}] .

Here the inhomogeneous part is represented as a sum of two vectors. Since the system is linear, we can apply the principle of superposition to find the solution.

First consider the homogeneous system:

X^{'} (t) = K X (t) .

Determine the eigenvalues of the matrix $K :$

det (K - λ I) = 0, \Rightarrow λ^{2} + (k_{1} + 2 k_{2} + k_{3} + k_{4}) λ + k_{1} k_{2} + k_{2} k_{3} + k_{3} k_{4} + k_{1} k_{4} + k_{2} k_{4} = 0.

Note that the coefficients in the resulting auxiliary equation, which is the quadratic equation, are always positive:

k_{1} + 2 k_{2} + k_{3} + k_{4} > 0, k_{1} k_{2} + k_{2} k_{3} + k_{3} k_{4} + k_{1} k_{4} + k_{2} k_{4} > 0.

Therefore, according to the Routh-Hurwitz criterion, the solution of the homogeneous system is asymptotically stable. In this case, it is the zero solution, since it is exactly the position of equilibrium of the homogeneous system.

Find the discriminant of the auxiliary equation:

D = b^{2} - 4 a c = {(k_{1} + 2 k_{2} + k_{3} + k_{4})}^{2} - 4 (k_{1} k_{2} + k_{2} k_{3} + k_{3} k_{4} + k_{1} k_{4} + k_{2} k_{4}) = k_{1}^{2} + 4 k_{2}^{2} + k_{3}^{2} + k_{4}^{2} + 4 k_{1} k_{2} + 2 k_{1} k_{3} + 2 k_{1} k_{4} + 4 k_{2} k_{3} + 4 k_{2} k_{4} + 2 k_{3} k_{4} - 4 k_{1} k_{2} - 4 k_{2} k_{3} - 4 k_{3} k_{4} - 4 k_{1} k_{4} - 4 k_{2} k_{4} = k_{1}^{2} + 4 k_{2}^{2} + k_{3}^{2} + k_{4}^{2} + 2 k_{1} k_{3} - 2 k_{1} k_{4} - 2 k_{3} k_{4} = {(k_{1} + k_{3})}^{2} + 4 k_{2}^{2} + k_{4}^{2} - 2 k_{1} k_{4} - 2 k_{3} k_{4} = {(k_{1} + k_{3})}^{2} - 2 k_{4} (k_{1} + k_{3}) + k_{4}^{2} + 4 k_{2}^{2} = {(k_{1} + k_{3} - k_{4})}^{2} + 4 k_{2}^{2} > 0.

Thus, the discriminant of the auxiliary equation is always positive. This means that the point of equilibrium is always a "node". As the stability of the system has been already proven above, this point is always a stable node, i.e. the eigenvalues $λ_{1}, λ_{2}$ are distinct and have real negative values. They are defined by the following equation:

λ_{1, 2} = \frac{1}{2} {- (k_{1} + 2 k_{2} + k_{3} + k_{4}) \pm {[{(k_{1} + k_{3} - k_{4})}^{2} + 4 k_{2}^{2}]}^{\frac{1}{2}}} .

Next, to avoid cumbersome expressions, we consider particular values of the thermal conductivity coefficients $k_{i} .$ Take, for example, the following values:

k_{1} = \frac{1}{10}, k_{2} = \frac{1}{5}, k_{3} = \frac{2}{5}, k_{4} = \frac{1}{2} .

The higher the value of $k_{i},$ the worse the insulation and the higher the thermal conductivity of the corresponding surface. With this choice of the coefficients $k_{i}$ the discriminant $D$ of the characteristic equation will be equal to

D = {(k_{1} + k_{3} - k_{4})}^{2} + 4 k_{2}^{2} = {(\frac{1}{10} + \frac{2}{5} - \frac{1}{2})}^{2} + 4 \cdot {(\frac{1}{5})}^{2} = \frac{4}{25} .

Then the eigenvalues $λ_{1, 2}$ have the following values:

λ_{1, 2} = \frac{- (k_{1} + 2 k_{2} + k_{3} + k_{4}) \pm \sqrt{D}}{2} = - \frac{9}{10}, - \frac{1}{2} .

Now we compute the eigenvectors associated with the eigenvalues $λ_{1}, λ_{2} .$ Let the eigenvalue $λ_{1} = - \frac{9}{10}$ be related to the vector $V_{1} = {(V_{11}, V_{21})}^{T},$ where the superscript $^{T}$ denotes transposition. We obtain:

(K - λ_{1} I) V_{1} = 0, \Rightarrow [\begin{matrix} - k_{1} - k_{2} - k_{3} - λ_{1} & k_{2} \\ k_{2} & - k_{2} - k_{4} - λ_{1} \end{matrix}] [\begin{matrix} V_{11} \\ V_{21} \end{matrix}] = 0, \Rightarrow [\begin{matrix} \frac{1}{5} & \frac{1}{5} \\ \frac{1}{5} & \frac{1}{5} \end{matrix}] [\begin{matrix} V_{11} \\ V_{21} \end{matrix}] = 0, \Rightarrow \frac{1}{5} V_{11} + \frac{1}{5} V_{21} = 0, \Rightarrow V_{11} = - V_{21}, \Rightarrow V_{1} = [\begin{array}{r} - 1 \\ 1 \end{array}] .

Similarly, we find the second eigenvector $V_{2} = {(V_{12}, V_{22})}^{T}$ associated with the eigenvalue $λ_{2} = - \frac{1}{2} :$

(K - λ_{2} I) V_{2} = 0, \Rightarrow [\begin{matrix} - \frac{1}{5} & \frac{1}{5} \\ \frac{1}{5} & - \frac{1}{5} \end{matrix}] [\begin{matrix} V_{12} \\ V_{22} \end{matrix}] = 0, \Rightarrow - \frac{1}{5} V_{12} + \frac{1}{5} V_{22} = 0, \Rightarrow V_{12} = V_{22}, \Rightarrow V_{2} = [\begin{array}{r} 1 \\ 1 \end{array}] .

The general solution of the homogeneous system with the given coefficients $k_{i}$ can be written as

X_{0} (t) = [\begin{matrix} x_{0} (t) \\ y_{0} (t) \end{matrix}] = C_{1} e^{- \frac{9}{10} t} [\begin{matrix} - 1 \\ 1 \end{matrix}] + C_{2} e^{- \frac{1}{2} t} [\begin{matrix} 1 \\ 1 \end{matrix}],

where the constants $C_{1},$ $C_{2}$ depend on the initial conditions.

Now we construct a solution of the nonhomogeneous system. First, consider the nonhomogeneous part $F_{1},$ which describes the parameters of the external environment:

F_{1} (t) = [\begin{matrix} k_{1} T_{g} + k_{3} T_{e} (t) \\ k_{4} T_{e} (t) \end{matrix}] = [\begin{matrix} \frac{1}{10} T_{g} + \frac{2}{5} T_{e} (t) \\ \frac{1}{2} T_{e} (t) \end{matrix}] = \frac{1}{10} [\begin{matrix} T_{g} + 4 T_{e} (t) \\ 5 T_{e} (t) \end{matrix}] .

Suppose, for example, the temperature of the ground (we will measure it in degrees Celsius) is $T_{g} = 10^{\circ} C .$ Assume also that the ambient temperature changes with daily periodicity by the law of

T_{e} (t) = 10 + 10 \sin (\frac{2 π}{24} t) = 10 + 10 \sin (\frac{π}{12} t),

that is, in the range from $0^{\circ} C$ to $20^{\circ} C .$ The time $t$ and the period of oscillations are expressed in hours. Then we have:

F_{1} (t) = \frac{1}{10} [\begin{matrix} 10 + 4 (10 + 10 \sin (\frac{π}{12} t)) \\ 5 (10 + 10 \sin (\frac{π}{12} t)) \end{matrix}] = [\begin{matrix} 5 + 4 \sin (\frac{π}{12} t) \\ 5 + 5 \sin (\frac{π}{12} t) \end{matrix}] .

We seek a particular solution $X_{1} (t)$ for this nonhomogeneous part in the form of vector quasipolynomial

X_{1} (t) = [\begin{matrix} x_{1} (t) \\ y_{1} (t) \end{matrix}] = [\begin{matrix} D_{1} + A_{1} \sin (\frac{π}{12} t) + B_{1} \cos (\frac{π}{12} t) \\ D_{2} + A_{2} \sin (\frac{π}{12} t) + B_{2} \cos (\frac{π}{12} t) \end{matrix}] .

Substitute the trial solution $X_{1} (t)$ into the system of equations with the inhomogeneous part $F_{1} (t) :$

{X^{'}}_{1} (t) = K X_{1} (t) + F_{1} (t),

\Rightarrow [\begin{matrix} x_{1}^{'} (t) \\ y_{1}^{'} (t) \end{matrix}] = [\begin{matrix} - k_{1} - k_{2} - k_{3} & k_{2} \\ k_{2} & - k_{2} - k_{4} \end{matrix}] [\begin{matrix} x_{1} (t) \\ y_{1} (t) \end{matrix}] + [\begin{matrix} k_{1} T_{g} + k_{3} T_{e} (t) \\ k_{4} T_{e} (t) \end{matrix}],

\Rightarrow [\begin{matrix} x_{1}^{'} (t) \\ y_{1}^{'} (t) \end{matrix}] = [\begin{matrix} - \frac{7}{10} & \frac{2}{10} \\ \frac{2}{10} & - \frac{7}{10} \end{matrix}] [\begin{matrix} x_{1} (t) \\ y_{1} (t) \end{matrix}] + [\begin{matrix} 5 + 4 \sin (\frac{π}{12} t) \\ 5 + 5 \sin (\frac{π}{12} t) \end{matrix}],

As a result we have

A_{1} \frac{π}{12} \cos (\frac{π}{12} t) - B_{1} \frac{π}{12} \sin (\frac{π}{12} t) = - \frac{7}{10} D_{1} - \frac{7}{10} A_{1} \sin (\frac{π}{12} t) - \frac{7}{10} B_{1} \cos (\frac{π}{12} t) + \frac{2}{10} D_{2} + \frac{2}{10} A_{2} \sin (\frac{π}{12} t) + \frac{2}{10} B_{2} \cos (\frac{π}{12} t) + 5 + 4 \sin (\frac{π}{12} t),

A_{2} \frac{π}{12} \cos (\frac{π}{12} t) - B_{2} \frac{π}{12} \sin (\frac{π}{12} t) = \frac{2}{10} D_{1} + \frac{2}{10} A_{1} \sin (\frac{π}{12} t) + \frac{2}{10} B_{1} \cos (\frac{π}{12} t) - \frac{7}{10} D_{2} - \frac{7}{10} A_{2} \sin (\frac{π}{12} t) - \frac{7}{10} B_{2} \cos (\frac{π}{12} t) + 5 + 5 \sin (\frac{π}{12} t) .

By equating the coefficients of like terms, we obtain the following algebraic system:

{\begin{cases} A_{1} \frac{π}{12} = - \frac{7}{10} B_{1} + \frac{2}{10} B_{2} \\ - B_{1} \frac{π}{12} = - \frac{7}{10} A_{1} + \frac{2}{10} A_{2} + 4 \\ 0 = - \frac{7}{10} D_{1} + \frac{2}{10} D_{2} + 5 \\ A_{2} \frac{π}{12} = \frac{2}{10} B_{1} - \frac{7}{10} B_{2} \\ - B_{2} \frac{π}{12} = \frac{2}{10} A_{1} - \frac{7}{10} A_{2} + 5 \\ 0 = \frac{2}{10} D_{1} - \frac{7}{10} D_{2} + 5 \end{cases} .

Solving this system, we find the coefficients $A_{1},$ $B_{1},$ $D_{1},$ $A_{2},$ $B_{2},$ $D_{2} :$

A_{1} = 6.55, B_{1} = - 3.55, A_{2} = 7.58, B_{2} = - 3.85, D_{1} = D_{2} = 10.

Thus, the particular solution $X_{1} (t)$ can be written as

X_{1} (t) = [\begin{matrix} x_{1} (t) \\ y_{1} (t) \end{matrix}] = [\begin{matrix} 10 + 6.55 \sin (\frac{π t}{12}) - 3.55 \cos (\frac{π t}{12}) \\ 10 + 7.58 \sin (\frac{π t}{12}) - 3.85 \cos (\frac{π t}{12}) \end{matrix}] .

Similarly, we find a particular solution $X_{2} (t)$ corresponding to the inhomogeneous part $F_{2} (t) :$

F_{2} (t) = [\begin{matrix} f (t) \\ 0 \end{matrix}] .

The heat source $f (t)$ can be programmed to different modes of operation. We restrict ourselves to the simplest case, when the heat source operates on a constant level: $f (t) = f_{0} .$ Then we seek a particular solution $X_{2}$ as a constant vector:

X_{2} = [\begin{matrix} G \\ H \end{matrix}] .

After substitution of the solution into the system with the inhomogeneous part $F_{2}$ we obtain:

{X^{'}}_{2} = K X_{2} + F_{2}, \Rightarrow 0 = [\begin{matrix} - \frac{7}{10} & \frac{2}{10} \\ \frac{2}{10} & - \frac{7}{10} \end{matrix}] [\begin{matrix} G \\ H \end{matrix}] + [\begin{matrix} f_{0} \\ 0 \end{matrix}], \Rightarrow {\begin{matrix} - \frac{7}{10} G + \frac{2}{10} H + f_{0} = 0 \\ \frac{2}{10} G - \frac{7}{10} H = 0 \end{matrix}, \Rightarrow {\begin{matrix} - 7 G + 2 H + 10 f_{0} = 0 \\ 2 G - 7 H = 0 \end{matrix}, \Rightarrow H = \frac{4}{9} f_{0}, \Rightarrow G = \frac{14}{9} f_{0} .

Thus, the particular solution $X_{2}$ is represented by the following constant vector:

X_{2} = [\begin{matrix} x_{2} \\ y_{2} \end{matrix}] = [\begin{matrix} G \\ H \end{matrix}] = [\begin{matrix} \frac{14}{9} f_{0} \\ \frac{4}{9} f_{0} \end{matrix}] .

According to the principle of superposition, a particular solution of the linear system, the inhomogeneous part of which is the sum of the vectors $F_{1} (t) + F_{2},$ is expressed as the sum of the corresponding particular solutions $X_{1} (t) + X_{2} .$ As a result, the general solution can be written as

X (t) = [\begin{matrix} x (t) \\ y (t) \end{matrix}] = X_{0} (t) + X_{1} (t) + X_{2} = C_{1} e^{- \frac{9}{10} t} [\begin{array}{r} - 1 \\ 1 \end{array}] + C_{2} e^{- \frac{1}{2} t} [\begin{matrix} 1 \\ 1 \end{matrix}] + [\begin{matrix} 10 + 6.55 \sin (\frac{π t}{12}) - 3.55 \cos (\frac{π t}{12}) \\ 10 + 7.58 \sin (\frac{π t}{12}) - 3.85 \cos (\frac{π t}{12}) \end{matrix}] + [\begin{matrix} \frac{14}{9} f_{0} \\ \frac{4}{9} f_{0} \end{matrix}] .

Find the constants $C_{1},$ $C_{2}$ from the initial conditions. Suppose, for example, $x (0) = 10^{\circ} C,$ $y (0) = 5^{\circ} C .$ Consequently,

{\begin{cases} x (0) = - C_{1} + C_{2} + 10 - 3.55 + \frac{14}{9} f_{0} = 10 \\ y (0) = C_{1} + C_{2} + 10 - 3.85 + \frac{4}{9} f_{0} = 5 \end{cases}, \Rightarrow {\begin{array}{l} - C_{1} + C_{2} = 3.55 - 1.56 f_{0} \\ C_{1} + C_{2} = - 1.55 - 0.44 f_{0} \end{array}, \Rightarrow {\begin{array}{l} C_{1} = - 2.35 + 0.56 f_{0} \\ C_{2} = 1.20 - f_{0} \end{array} .

The final general solution given the initial conditions and for an arbitrary constant heat source is written in the form:

X (t) = [\begin{matrix} x (t) \\ y (t) \end{matrix}] = (- 2.35 + 0.56 f_{0}) e^{- \frac{9}{10} t} [\begin{array}{r} - 1 \\ 1 \end{array}] + (1.20 - f_{0}) e^{- \frac{1}{2} t} [\begin{matrix} 1 \\ 1 \end{matrix}] + [\begin{matrix} 10 + 6.55 \sin (\frac{π t}{12}) - 3.55 \cos (\frac{π t}{12}) \\ 10 + 7.58 \sin (\frac{π t}{12}) - 3.85 \cos (\frac{π t}{12}) \end{matrix}] + [\begin{matrix} 1.56 \\ 0.44 \end{matrix}] f_{0} .

Here, the first two terms are damped and describe the transition process. The third term is related to heat exchange with the environment. Recall that in this problem the following law of the change of the ambient temperature has been chosen:

T_{e} (t) = 10 + 10 \sin (\frac{π}{12} t) .

Finally, the last term in the formula describes the action of the constant heat source $f_{0} .$ The parameter $f_{0}$ is proportional to the source power. The heat source working for an hour provides an increase in temperature (on the ground floor) by $f_{0}$ degrees.

Consider the graphs of temperature change. Figure $2$ shows the temperature curves on the lower and upper floor of the house (respectively, the functions $x (t)$ and $y (t)$ ) when the heater is turned off: $f_{0} = 0.$ It can be seen that an increase of the outside temperature $T_{e} (t)$ causes the temperature rise inside the house with a delay. The inertia of heat transfer also causes the temperature on the ground floor to be changed in the range from ${2.5}^{\circ} C$ to ${17.5}^{\circ} C,$ and on the top floor − from ${1.5}^{\circ} C$ to ${18.5}^{\circ} C,$ while the ambient temperature ranges from $0^{\circ} C$ to $20^{\circ} C .$

When the heat source $f_{0}$ is turned on, of course, it is getting warmer in the house. Figure $3$ shows the graphs of $x (t),$ $y (t)$ for the case $f_{0} = 5 (\frac{deg}{hr}) .$ Now, the temperature on the ground floor ranges from $10^{\circ} C$ to $25^{\circ} C,$ and on the upper floor − from $5^{\circ} C$ to ${20.5}^{\circ} C .$

As the power of the heat source increases, the ground floor is heated much better than the upper floor (Figure $4$ ).

Exploring the solutions $x (t),$ $y (t)$ under different operating conditions of the heat source, you can choose the most comfortable and cost effective heating mode.