8  Thermal Dynamics of Buildings: Part II

8.1 Lecture Overview

NoteLearning Objectives

By the end of this module, students will be able to:

• Derive thermal network models from first principles:
  • Describe the elements of the Fourier heat equation (the 1-D Fourier transient heat conduction PDE and its steady-state version), including its solution as a sum of steady state plus an infinite sum of exponentials
  • Recognize the link between the Fourier heat equation PDE and the Thermal Network Model, understanding it as a finite difference approach using physical construction-based slicing/lumping that leads to a set of N simultaneous first-order differential equations to be solved
  • Understand the link between Kirchhoff’s voltage/current laws and the heat balance equations we can work out in thermal network models
  • Recognize that the linearity assumption means with respect to the behavior of network elements (independent of temperature and time), and other implications
  • Directly represent radiation heat transfer interactions, lumped building mass components, solar radiation, heat extraction by HVAC, infiltration heat paths
• Solve 1R1C and 2R2C networks to find a building’s time constants
  • Understand the importance of time constants even if they are just approximations
• Write down the linear state-space representation of an arbitrary thermal network model

NoteTopics Covered

• The Fourier heat equation: from steady-state to transient
• Analytical solution of the heat equation
• Finite difference discretization and lumped parameter models
• From PDEs to thermal network models (RC circuits)
• Kirchhoff’s laws and heat balance equations
• Linearity assumptions and their implications
• Representing heat transfer elements in networks
• Solving 1R1C networks: the single time constant model
• Solving 2R2C networks: multiple time constants
• Physical interpretation of time constants
• State-space representation of thermal networks

TipProject Milestones

This module provides the mathematical foundation for building thermal models that can be:

• Used in model predictive control (MPC) strategies
• Identified from data using system identification techniques
• Simulated to predict building response to weather and HVAC inputs
• Analyzed for controllability and observability

Understanding these models is essential for designing intelligent HVAC controllers in your final project.

8.2 From Steady-State to Dynamic: The Need for Time

In Lecture 5, we developed a solid understanding of steady-state heat transfer. We learned to calculate heat loss rates, R-values, and U-values—but we largely ignored time. The questions we posed at the end of that lecture highlighted this gap:

• How long does it take for a building to respond to temperature changes?
• When should we start heating before occupants arrive?
• Why do massive buildings “coast” through cold snaps better than lightweight ones?

To answer these questions, we need to understand how temperatures evolve over time. This requires moving from algebraic equations to differential equations.

8.2.1 A Motivating Example: The Office Room Revisited

Let’s revisit a concrete example to frame the questions we’ll answer in this lecture. Consider a small office room with dimensions 4m \(\times\) 4m \(\times\) 2.7m. The room has one exterior wall (4m \(\times\) 2.7m = 10.8 m\(^2\)) containing a 2 m\(^2\) double-pane window. The other walls, floor, and ceiling are adjacent to conditioned spaces (assume adiabatic for this analysis).

Given:

• Exterior wall U-value: 0.35 W/m\(^2\)\(\cdot\)K
• Window U-value: 2.0 W/m\(^2\)\(\cdot\)K
• Air density: 1.2 kg/m\(^3\)
• Air specific heat: 1,005 J/kg\(\cdot\)K
• Indoor temperature: 20\(^\circ\)C
• Outdoor temperature: 0\(^\circ\)C

From Lecture 5, we can calculate the steady-state heat loss:

• Wall conductance: \(U_{wall} \times A_{wall} = 0.35 \times 8.8 = 3.08\) W/K
• Window conductance: \(U_{window} \times A_{window} = 2.0 \times 2.0 = 4.00\) W/K
• Total UA: 7.08 W/K
• Heat loss: \(\dot{Q} = UA \times \Delta T = 7.08 \times 20 = 142\) W

But now let’s ask dynamic questions:

1. What is the thermal capacitance of the air in this room? \[C_{air} = \rho \cdot V \cdot c = 1.2 \times 43.2 \times 1005 = 52{,}100 \text{ J/K} \approx 52 \text{ kJ/K}\]

2. What is the time constant of this room? \[\tau = R \times C = \frac{C}{UA} = \frac{52{,}100}{7.08} = 7{,}360 \text{ s} \approx 2 \text{ hours}\]

3. If the heating fails at midnight when \(T_{in} = 20^\circ\text{C}\) and \(T_{out} = 0^\circ\text{C}\), how does the temperature evolve? How long until it drops below 15\(^\circ\)C?

4. If this room had a 50mm thick exposed concrete floor (16 m\(^2\), \(\rho = 2300\) kg/m\(^3\), \(c = 880\) J/kg\(\cdot\)K), the thermal capacitance would jump to over 4,000 kJ/K, and the time constant would become ~73 hours. Why does this matter for building control?

By the end of this lecture, you’ll be able to answer all of these questions and understand the mathematical framework behind them. Notice how the window, despite being only 19% of the wall area, contributes 57% of the heat loss—a reminder that thermal networks can reveal non-obvious insights about building performance.

8.3 The Fourier Heat Equation

8.3.1 The 1-D Transient Heat Conduction PDE

In Lecture 5, we derived Fourier’s law for steady-state conduction:

\[\dot{Q} = -kA\frac{dT}{dx}\]

But this assumes temperatures don’t change with time. To understand dynamic behavior, we need to derive the transient heat equation from an energy balance.

Consider a thin slab of material with thickness \(\Delta x\), area \(A\), density \(\rho\), and specific heat \(c\). The rate of change of energy stored in this slab equals the net heat flow into it:

\[\rho c A \Delta x \frac{\partial T}{\partial t} = \dot{Q}_{in} - \dot{Q}_{out}\]

Using Fourier’s law for the heat flows at positions \(x\) and \(x + \Delta x\):

\[\rho c A \Delta x \frac{\partial T}{\partial t} = -kA\frac{\partial T}{\partial x}\bigg|_x - \left(-kA\frac{\partial T}{\partial x}\bigg|_{x+\Delta x}\right)\]

Dividing by \(A \Delta x\) and taking the limit as \(\Delta x \to 0\):

\[\rho c \frac{\partial T}{\partial t} = k\frac{\partial^2 T}{\partial x^2}\]

TipConvince yourself: Algebraic steps to derive the heat equation

Let’s work through the algebra carefully to see how we get from the energy balance to the heat equation.

Starting equation: \[\rho c A \Delta x \frac{\partial T}{\partial t} = -kA\frac{\partial T}{\partial x}\bigg|_x - \left(-kA\frac{\partial T}{\partial x}\bigg|_{x+\Delta x}\right)\]

Step 1: Handle the double negative

The second term has \(-(-kA...)\), which becomes positive: \[\rho c A \Delta x \frac{\partial T}{\partial t} = -kA\frac{\partial T}{\partial x}\bigg|_x + kA\frac{\partial T}{\partial x}\bigg|_{x+\Delta x}\]

Step 2: Factor out \(kA\) and rearrange

\[\rho c A \Delta x \frac{\partial T}{\partial t} = kA\left(\frac{\partial T}{\partial x}\bigg|_{x+\Delta x} - \frac{\partial T}{\partial x}\bigg|_x\right)\]

Note the order: we have (gradient at \(x+\Delta x\)) minus (gradient at \(x\)).

Step 3: Divide by \(A \Delta x\)

\[\rho c \frac{\partial T}{\partial t} = \frac{k}{\Delta x}\left(\frac{\partial T}{\partial x}\bigg|_{x+\Delta x} - \frac{\partial T}{\partial x}\bigg|_x\right)\]

Step 4: Take the limit as \(\Delta x \to 0\)

The term in brackets divided by \(\Delta x\) is exactly the definition of the derivative of \(\frac{\partial T}{\partial x}\) with respect to \(x\):

\[\lim_{\Delta x \to 0}\frac{1}{\Delta x}\left(\frac{\partial T}{\partial x}\bigg|_{x+\Delta x} - \frac{\partial T}{\partial x}\bigg|_x\right) = \frac{\partial^2 T}{\partial x^2}\]

Final result: \[\rho c \frac{\partial T}{\partial t} = k\frac{\partial^2 T}{\partial x^2}\]

Or, in its standard form:

\[\boxed{\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}}\]

where \(\alpha = k/(\rho c)\) is the thermal diffusivity (units: m\(^2\)/s).

Physical interpretation:

• The left side, \(\partial T/\partial t\), is the rate of temperature change at a point
• The right side, \(\alpha \partial^2 T/\partial x^2\), involves the curvature of the temperature profile
• Where the temperature profile is curved (concave up), heat flows in from both sides, causing temperature to rise
• Where the profile is straight (no curvature), there’s no net accumulation of heat

TipConvince yourself: What does “curvature” mean physically?

The second derivative \(\partial^2 T/\partial x^2\) measures how much the temperature gradient changes with position. Consider these three scenarios through a wall:

1. Hot spot in the middle (temperature profile curves downward, concave down):

T |     *  ← hot spot
  |    / \
  |   /   \
  |  /     \
  |---------→ x

Heat flows away from the hot spot toward both cooler sides. The gradient \(\partial T/\partial x\) is positive (upward slope) on the left and negative (downward slope) on the right. As we move through the hot spot, the gradient decreases (becomes less positive, then negative), so \(\partial^2 T/\partial x^2 < 0\). Result: temperature at the hot spot decreases over time.

2. Cold spot in the middle (temperature profile curves upward, concave up):

T |  \     /
  |   \   /
  |    \ /
  |     * ← cold spot
  |---------→ x

Heat flows toward the cold spot from both warmer sides. The gradient is negative (downward slope) on the left and positive (upward slope) on the right. As we move through the cold spot, the gradient increases, so \(\partial^2 T/\partial x^2 > 0\). Result: temperature at the cold spot increases over time.

3. Linear profile (no curvature):

T |  \
  |   \
  |    \
  |     \
  |---------→ x

The gradient is constant throughout. No matter where we are, heat flows through at a steady rate—what comes in equals what goes out. \(\partial^2 T/\partial x^2 = 0\). Result: temperature doesn’t change with time (steady state).

Key insight: The heat equation says temperature changes fastest where curvature is greatest. Heat “diffuses” from convex regions (hot spots) toward concave regions (cold spots), gradually smoothing out the temperature profile until it becomes linear (steady state).

8.3.2 Steady-State as a Special Case

When the system reaches steady state, temperatures stop changing: \(\partial T/\partial t = 0\). The heat equation reduces to:

\[\frac{d^2 T}{dx^2} = 0\]

This equation says the temperature profile has zero curvature—it must be a straight line! This is exactly what we assumed in Lecture 5 when we wrote:

\[\dot{Q} = kA\frac{T_1 - T_2}{\Delta x}\]

The linear temperature profile through a wall is the steady-state solution to the Fourier heat equation.

TipConvince yourself

Let’s show explicitly how the zero-curvature condition leads to the Lecture 5 formula.

Starting from the steady-state heat equation:

\[\frac{d^2 T}{dx^2} = 0\]

Solving this ordinary differential equation:

First integration gives: \[\frac{dT}{dx} = C_1\]

where \(C_1\) is a constant. Second integration gives: \[T(x) = C_1 x + C_2\]

This confirms the temperature profile is linear in space.

Applying boundary conditions:

Consider a wall extending from \(x=0\) to \(x=\Delta x\), with temperatures \(T(0) = T_1\) (hot side) and \(T(\Delta x) = T_2\) (cold side).

From \(T(0) = T_1\): \[C_2 = T_1\]

From \(T(\Delta x) = T_2\): \[T_2 = C_1 \Delta x + T_1\] \[C_1 = \frac{T_2 - T_1}{\Delta x}\]

The key insight:

Recall from our first integration that \(\frac{dT}{dx} = C_1\) everywhere in the wall (the gradient is constant). We’ve now determined from the boundary conditions that \(C_1 = \frac{T_2 - T_1}{\Delta x}\). Therefore, by substitution:

\[\frac{dT}{dx} = \frac{T_2 - T_1}{\Delta x} \quad \text{(constant throughout the wall)}\]

Applying Fourier’s law:

Since the gradient is constant everywhere in the wall, we can apply Fourier’s law: \[\dot{Q} = -kA\frac{dT}{dx} = -kA\frac{T_2 - T_1}{\Delta x} = kA\frac{T_1 - T_2}{\Delta x}\]

This is exactly the formula we used in Lecture 5! It’s not an approximation or assumption—it’s the exact steady-state solution to the Fourier heat equation. The linear temperature profile with constant gradient emerges naturally when \(\partial T/\partial t = 0\).

8.3.3 Solution Structure: Steady-State Plus Transients

The general solution to the heat equation (with appropriate boundary conditions) can be found using a technique called separation of variables. This method splits the PDE into two ODEs—one for time, one for space. The time component leads to exponential decay, while the spatial component gives rise to eigenfunctions. The result is a beautiful structure:

\[T(x,t) = T_{ss}(x) + \sum_{n=1}^{\infty} A_n \phi_n(x) e^{-t/\tau_n}\]

where:

• \(T_{ss}(x)\) is the steady-state solution—the linear profile we calculated in Lecture 5
• The summation represents transient modes, each with its own spatial shape \(\phi_n(x)\) and time constant \(\tau_n\)
• The coefficients \(A_n\) depend on the initial conditions

Key insight: As \(t \to \infty\), all the exponential terms decay to zero, leaving only the steady-state solution. The transient terms describe how the system approaches steady state, not where it ends up.

TipExplicit form with sinusoidal eigenfunctions

The abstract form \(T(x,t) = T_{ss}(x) + \sum_{n=1}^{\infty} A_n \phi_n(x) e^{-t/\tau_n}\) hides the spatial details. Let’s connect it to the explicit solution with sinusoidal eigenfunctions that emerges from separation of variables.

The most general solution to the 1-D heat equation can be written as:

\[T(x,t) = A_0 + B_0 x + \sum_{n=1}^{\infty} \big[ A_n \cos(\lambda_n x) + B_n \sin(\lambda_n x) \big] e^{-\lambda_n^2 \alpha t}\]

where \(\lambda_n\) are eigenvalues determined by the boundary conditions.

Connecting the two forms:

1. Steady-state part: \(T_{ss}(x) \leftrightarrow A_0 + B_0 x\)
   • For fixed boundary temperatures, this is the linear profile we derived earlier
   • \(B_0 = (T_2 - T_1)/L\) gives the constant gradient
   • \(A_0 = T_1\) sets the reference temperature
2. Spatial eigenfunctions: \(\phi_n(x) \leftrightarrow A_n \cos(\lambda_n x) + B_n \sin(\lambda_n x)\)
   • These are the sinusoidal “shapes” that the temperature profile can take
   • For a wall with fixed temperatures at \(x=0\) and \(x=L\), boundary conditions force \(A_n = 0\) (no cosine terms) and \(\lambda_n = n\pi/L\)
   • This gives \(\phi_n(x) = \sin(n\pi x/L)\)
3. Time constants: \(\tau_n = \frac{1}{\lambda_n^2 \alpha}\)
   • The exponential \(e^{-t/\tau_n}\) is exactly the same as \(e^{-\lambda_n^2 \alpha t}\)
   • For \(\lambda_n = n\pi/L\): \(\tau_n = \frac{L^2}{n^2 \pi^2 \alpha}\)
   • Higher modes (\(n=2, 3, \ldots\)) decay as \(1/n^2\), so they vanish quickly

Physical meaning of the eigenfunctions:

Each \(\sin(n\pi x/L)\) represents a spatial “mode” with \(n\) half-wavelengths across the wall thickness:

• \(n=1\): One smooth arc from \(x=0\) to \(x=L\) (fundamental mode, slowest decay)
• \(n=2\): Two half-waves (decays 4× faster)
• \(n=3\): Three half-waves (decays 9× faster)

The specific mix of these modes (the coefficients \(A_n\) or \(B_n\)) depends on the initial temperature distribution \(T(x,0)\), determined by Fourier analysis.

The time constants \(\tau_n\) tell us how quickly each mode decays. From the eigenvalue analysis (detailed in the derivation below), for a wall of thickness \(L\):

\[\tau_n \propto \frac{L^2}{\alpha n^2}\]

The first mode (\(n=1\)) decays slowest and dominates the long-term response. Higher modes decay much faster (proportional to \(1/n^2\)).

TipConvince yourself: Separation of variables derivation

Let’s derive this solution structure from first principles using the method of separation of variables.

Step 1: Assume a separable solution

We look for solutions of the form: \[T(x,t) = X(x) \cdot \Theta(t)\]

where \(X(x)\) depends only on position and \(\Theta(t)\) depends only on time.

Step 2: Substitute into the heat equation

Starting with \(\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}\) and substituting our assumed form:

\[X(x) \frac{d\Theta}{dt} = \alpha \Theta(t) \frac{d^2X}{dx^2}\]

Step 3: Separate the variables

Divide both sides by \(X(x)\Theta(t)\):

\[\frac{1}{\Theta}\frac{d\Theta}{dt} = \frac{\alpha}{X}\frac{d^2X}{dx^2}\]

The left side depends only on \(t\), the right side only on \(x\). For this to be true for all \(x\) and \(t\), both sides must equal the same constant. We call this constant \(-\lambda\) (the negative sign is chosen for convenience):

\[\frac{1}{\Theta}\frac{d\Theta}{dt} = -\lambda \quad \text{and} \quad \frac{\alpha}{X}\frac{d^2X}{dx^2} = -\lambda\]

Step 4: Solve the time equation

From \(\frac{d\Theta}{dt} = -\lambda \Theta\), we get: \[\Theta(t) = e^{-\lambda t}\]

If we define \(\tau = 1/\lambda\), this becomes \(\Theta(t) = e^{-t/\tau}\)—our exponential decay!

Step 5: Solve the spatial equation

The spatial equation becomes: \[\frac{d^2X}{dx^2} = -\frac{\lambda}{\alpha}X\]

This is an eigenvalue problem. The solutions are sinusoidal functions whose specific form depends on the boundary conditions. For a wall with fixed boundary temperatures, we get: \[X_n(x) = \phi_n(x) = \sin\left(\frac{n\pi x}{L}\right)\]

with corresponding eigenvalues: \[\lambda_n = \frac{n^2 \pi^2 \alpha}{L^2}\]

This gives time constants: \[\tau_n = \frac{1}{\lambda_n} = \frac{L^2}{n^2 \pi^2 \alpha}\]

Step 6: Superposition principle

Because the heat equation is linear, we can sum all possible solutions: \[T(x,t) = \sum_{n=1}^{\infty} A_n \phi_n(x) e^{-t/\tau_n}\]

Step 7: Add the steady-state solution

The complete solution must also satisfy the boundary conditions at all times. The full solution is: \[T(x,t) = T_{ss}(x) + \sum_{n=1}^{\infty} A_n \phi_n(x) e^{-t/\tau_n}\]

where \(T_{ss}(x)\) is the steady-state profile (linear, as we showed earlier), and the coefficients \(A_n\) are determined by matching the initial temperature distribution at \(t=0\).

Physical interpretation:

• Each mode \(n\) has its own spatial pattern \(\phi_n(x)\) and decay rate \(1/\tau_n\)
• Higher modes (\(n > 1\)) decay much faster than the fundamental mode (\(n=1\))
• After a few multiples of \(\tau_1\), only the steady-state solution remains
• The separation constant \(\lambda\) emerges as a physical eigenvalue that determines both the spatial pattern and the time constant

8.3.4 Thermal Diffusivity and Its Meaning

Thermal diffusivity \(\alpha = k/(\rho c)\) tells us how quickly temperature disturbances propagate through a material. It combines:

• \(k\): how easily heat conducts (higher \(k\) \(\rightarrow\) faster propagation)
• \(\rho c\): how much heat the material stores per unit volume (higher \(\rho c\) \(\rightarrow\) slower response)

Materials with high diffusivity respond quickly to surface temperature changes; materials with low diffusivity respond slowly.

Material	\(k\) (W/m\(\cdot\)K)	\(\rho\) (kg/m\(^3\))	\(c\) (J/kg\(\cdot\)K)	\(\alpha\) (mm\(^2\)/s)
Copper	385	8,940	385	112
Aluminum	205	2,700	900	84
Steel	50	7,850	460	14
Stone (granite)	2.8	2,600	820	1.3
Concrete	1.4	2,300	880	0.69
Brick	0.72	1,700	840	0.50
Glass	0.9	2,500	840	0.43
Gypsum board	0.16	800	1,000	0.20
Wood (softwood)	0.12	500	1,200	0.20
Water	0.6	1,000	4,186	0.14

Sources: Engineering ToolBox, Engineers Edge

Notice the huge range: copper’s diffusivity is nearly 1000× that of water! Metals respond almost instantaneously to surface heating; water and masonry materials respond very slowly.

8.3.5 Thermal Penetration Depth

A useful concept for building intuition is the thermal penetration depth—the distance a temperature disturbance travels into a material in a given time.

Recall from our eigenvalue analysis that the time constant scales as \(\tau \sim L^2/\alpha\), where \(L\) is a characteristic length. Rearranging this relationship, we can express the length scale that responds over a time \(t\) as \(L \sim \sqrt{\alpha t}\). This gives us the penetration depth.

For a transient event of duration \(t\):

\[\delta \approx \sqrt{\alpha t}\]

For periodic heating (like the daily outdoor temperature cycle with period \(T\)):

\[\delta \approx \sqrt{\frac{\alpha T}{\pi}}\]

This tells us how much of a wall’s mass actually “participates” in responding to temperature fluctuations.

TipLearn-by-Doing Activity: Thermal Penetration Depth

Using the formula \(\delta = \sqrt{\alpha t}\), calculate how far a temperature disturbance penetrates into different materials over various time scales.

Given thermal diffusivities:

• Concrete: \(\alpha = 0.69\) mm\(^2\)/s (\(6.9 \times 10^{-7}\) m\(^2\)/s)
• Brick: \(\alpha = 0.50\) mm\(^2\)/s
• Wood: \(\alpha = 0.20\) mm\(^2\)/s

Tasks:

1. Calculate the penetration depth for each material after:

   • 1 hour (3,600 s)
   • 1 day (86,400 s)
   • 1 week (604,800 s)

2. For daily outdoor temperature cycles, how deep does the temperature wave penetrate into a concrete wall? (Use \(\delta = \sqrt{\alpha T/\pi}\) with \(T = 86,400\) s)

3. A typical concrete wall is 200 mm thick. Based on your calculation, does the entire wall participate in daily temperature swings, or only the outer portion?

4. What does this imply about modeling walls for daily HVAC control versus seasonal energy calculations?

Answers:

Material	1 hour	1 day	1 week
Concrete	50 mm	244 mm	647 mm
Brick	43 mm	209 mm	552 mm
Wood	27 mm	132 mm	348 mm

For daily cycles in concrete: \(\delta \approx 138\) mm. This means only the outer ~140 mm of a concrete wall responds significantly to daily temperature swings—the inner core stays at a nearly constant temperature. This is why we can often “lump” the thermal mass rather than modeling the full distributed system.

8.4 From PDEs to Lumped Parameter Models

8.4.1 Why Lump? The Practical Challenge

The Fourier heat equation is a partial differential equation (PDE)—temperature varies continuously in both space and time. This is a distributed parameter system with infinitely many degrees of freedom.

For practical purposes—simulation, control design, parameter identification—we need models with a finite number of states. The solution is to discretize the spatial dimension, converting the PDE into a system of ordinary differential equations (ODEs).

The lumped parameter approach divides the continuous system into a finite number of nodes, each assumed to have a uniform temperature. The temperature gradient between nodes drives heat flow; the thermal mass at nodes determines how quickly temperatures change.

8.4.2 Finite Difference Discretization

Consider a wall divided into \(N\) nodes at positions \(x_1, x_2, \ldots, x_N\) with spacing \(\Delta x\). At each interior node \(i\), we approximate the second derivative using the central difference formula:

\[\frac{\partial^2 T}{\partial x^2}\bigg|_i \approx \frac{T_{i+1} - 2T_i + T_{i-1}}{(\Delta x)^2}\]

Substituting into the heat equation:

\[\frac{\partial T_i}{\partial t} = \alpha \frac{T_{i+1} - 2T_i + T_{i-1}}{(\Delta x)^2}\]

This is now an ordinary differential equation for \(T_i(t)\)—it depends only on time, not on continuous spatial derivatives. We have one such equation for each node, giving us \(N\) coupled first-order ODEs.

Rearranging:

\[\frac{dT_i}{dt} = \frac{\alpha}{(\Delta x)^2}(T_{i-1} - 2T_i + T_{i+1})\]

\[= \frac{\alpha}{(\Delta x)^2}(T_{i-1} - T_i) + \frac{\alpha}{(\Delta x)^2}(T_{i+1} - T_i)\]

Each term looks like heat flow driven by a temperature difference—divided by thermal resistance, into a thermal capacitance!

8.4.3 Physical Construction-Based Slicing

Rather than using uniform grid spacing, building scientists typically slice walls at material boundaries: the air film, drywall, insulation, sheathing, and exterior surface each become separate regions.

This “physical” discretization has several advantages:

1. Material properties change at boundaries: Conductivity, density, and specific heat are constant within each layer but jump at interfaces
2. Fewer nodes needed: A 3-node model of a wall (interior surface, insulation midpoint, exterior surface) often captures the essential dynamics
3. Physical interpretation: Each node corresponds to something tangible (the air, the wall’s interior mass, the exterior surface)

The standard approach assigns:

• One temperature to each material layer (or to each “lump” of a thick layer)
• Resistances connecting adjacent nodes based on layer thickness and conductivity
• Capacitances at nodes based on the mass and specific heat of each lump

8.4.4 The Emergence of RC Networks

Let’s see how the discretized equations map to RC circuits. Consider three adjacent nodes with temperatures \(T_{i-1}\), \(T_i\), and \(T_{i+1}\), connected by thermal resistances \(R_{i-1,i}\) and \(R_{i,i+1}\), with node \(i\) having thermal capacitance \(C_i\).

The heat balance at node \(i\) (Kirchhoff’s Current Law applied to heat flow):

\[C_i \frac{dT_i}{dt} = \frac{T_{i-1} - T_i}{R_{i-1,i}} + \frac{T_{i+1} - T_i}{R_{i,i+1}}\]

This is identical in form to the equation for the voltage at a node in an RC circuit! The heat equation PDE has been transformed into a thermal network model.

     T_{i-1}               T_i                T_{i+1}
        o----[R_{i-1,i}]----o----[R_{i,i+1}]----o
                            |
                           === C_i
                            |
                           GND

For a uniform material with spacing \(\Delta x\):

• \(R = \frac{\Delta x}{kA}\) (resistance of each segment)
• \(C = \rho c A \Delta x\) (capacitance of each node’s mass)

And the product \(RC = \frac{\rho c (\Delta x)^2}{k} = \frac{(\Delta x)^2}{\alpha}\), which is exactly the time scale that appears in our discretized equation.

The key insight: By discretizing the Fourier PDE, we’ve converted a distributed parameter system into a lumped parameter system that can be analyzed using circuit theory. All the tools of electrical network analysis—Kirchhoff’s laws, series/parallel combinations, nodal analysis—apply directly to thermal networks.

8.5 Thermal Network Models

8.5.1 The Electrical-Thermal Analogy Revisited

In Lecture 5, we introduced the analogy between thermal and electrical systems for steady-state (resistive) analysis. Now we extend it to include dynamic (capacitive) elements:

Electrical Domain	Thermal Domain	Units
Voltage \(V\)	Temperature \(T\)	V \(\leftrightarrow\) K (or \(^\circ\)C)
Current \(I\)	Heat flow rate \(\dot{Q}\)	A \(\leftrightarrow\) W
Charge \(q\)	Heat (energy) \(Q\)	C \(\leftrightarrow\) J
Resistance \(R_e\)	Thermal resistance \(R\)	\(\Omega\) \(\leftrightarrow\) K/W
Capacitance \(C_e\)	Thermal capacitance \(C\)	F \(\leftrightarrow\) J/K
Ohm’s Law: \(I = \frac{V_1 - V_2}{R_e}\)	Fourier’s Law: \(\dot{Q} = \frac{T_1 - T_2}{R}\)
Capacitor: \(I = C_e \frac{dV}{dt}\)	Thermal mass: \(\dot{Q} = C \frac{dT}{dt}\)
Time constant: \(\tau = R_e C_e\)	Time constant: \(\tau = RC\)	s

The analogy is remarkably complete. Every concept from electrical circuit analysis has a thermal counterpart:

• Series resistances add: \(R_{total} = R_1 + R_2\) (layers of a wall)
• Parallel resistances combine: \(\frac{1}{R_{total}} = \frac{1}{R_1} + \frac{1}{R_2}\) (parallel heat paths)
• RC time constants characterize dynamic response
• Kirchhoff’s laws enforce conservation

8.5.2 Kirchhoff’s Laws for Heat Balance

The two fundamental laws of circuit analysis translate directly to thermal networks:

8.5.2.1 Kirchhoff’s Current Law (KCL) \(\rightarrow\) Conservation of Energy

The sum of heat flows into any node must equal zero (in steady state) or equal the rate of heat storage (in transient).

For a node with temperature \(T_i\), capacitance \(C_i\), connected to neighboring nodes and receiving heat input \(\dot{Q}_i\):

\[C_i \frac{dT_i}{dt} = \sum_j \frac{T_j - T_i}{R_{ij}} + \dot{Q}_i\]

This is our fundamental nodal equation. It says:

• Rate of energy storage (left side) = net heat flow in (right side)
• Heat flows from higher to lower temperature (like current from higher to lower voltage)

8.5.2.2 Kirchhoff’s Voltage Law (KVL) \(\rightarrow\) Temperature Consistency

Around any closed loop, the sum of temperature differences must equal zero.

\[\sum_{\text{loop}} \Delta T = 0\]

This is automatically satisfied if we define temperatures consistently at each node. It’s less commonly used directly in thermal analysis but ensures our network is physically meaningful.

8.5.2.3 From KCL to a System of ODEs

Writing KCL at each capacitance node produces one differential equation per node. For a network with \(n\) capacitance nodes:

\[C_1 \frac{dT_1}{dt} = \text{(sum of heat flows into node 1)}\] \[C_2 \frac{dT_2}{dt} = \text{(sum of heat flows into node 2)}\] \[\vdots\] \[C_n \frac{dT_n}{dt} = \text{(sum of heat flows into node n)}\]

This system of \(n\) first-order ODEs completely describes the thermal dynamics of the building.

8.5.3 The Linearity Assumption

Our thermal network models are linear because we assume:

1. Resistances are constant: \(R\) doesn’t depend on temperature or time
2. Capacitances are constant: \(C\) doesn’t depend on temperature or time
3. Heat flows are proportional to temperature differences: \(\dot{Q} = \Delta T / R\)

When is linearity valid?

• Conduction: Generally linear. Thermal conductivity \(k\) varies weakly with temperature for most building materials.
• Convection: Approximately linear for forced convection; natural convection is weakly nonlinear but often linearized.
• Thermal mass: Linear as long as there are no phase changes (melting/freezing).

When does linearity break down?

• Radiation: Fundamentally nonlinear (\(\dot{Q} \propto T^4\)). However, we showed in Lecture 5 that for small temperature differences, radiation can be linearized as \(\dot{Q} \approx h_r A \Delta T\) with \(h_r = 4\varepsilon\sigma T_m^3\).
• Phase change materials: Highly nonlinear at the melting point.
• Variable convection: If airflow patterns change significantly with temperature.

Why linearity matters:

Linear systems have the superposition property: the response to multiple inputs equals the sum of responses to each input separately. This enables:

• Analytical solutions (eigenvalue methods)
• Efficient simulation
• Well-developed control theory (transfer functions, state-space methods)
• System identification from input-output data

For most building thermal analysis, the linearity assumption is excellent for temperature variations of $\(10-20\)^$C around typical operating points.

8.6 Representing Heat Transfer Elements

Now that we understand the general framework, let’s see how specific heat transfer mechanisms are represented in thermal networks.

8.6.1 Conduction Through Walls

From Lecture 5, we know that conduction through a layer of material is represented by a thermal resistance:

\[R_{cond} = \frac{\Delta x}{kA}\]

For a multi-layer wall, we can represent it as a chain of R-C-R-C elements:

T_out    R_se    T_s,o    R_1    T_1    R_2    T_2   R_si    T_in
  o------[===]-----o-----[===]----o----[===]----o----[===]----o
                   |             |             |
                  === C_1       === C_2       === C_in
                   |             |             |
                  GND           GND           GND

Where:

• \(R_{se}\), \(R_{si}\) are surface film resistances (convection + radiation)
• \(R_1\), \(R_2\) are conduction resistances of material layers
• \(C_1\), \(C_2\) are thermal capacitances of each layer

How many nodes are needed? The thermal penetration depth gives guidance:

• For fast dynamics (hourly control): Only the outer ~50-100 mm of massive materials participates. A 2-3 node model is usually sufficient.
• For slow dynamics (daily/seasonal): More of the mass participates. May need 3-5 nodes for thick concrete walls.
• For lightweight construction (wood frame with insulation): 1-2 nodes often suffice because there’s little thermal mass.

8.6.2 Convection at Surfaces

Surface heat transfer (combined convection and radiation to surroundings) is represented by a surface film resistance:

\[R_s = \frac{1}{h_c A}\]

where \(h_c\) is the combined surface heat transfer coefficient (convection + linearized radiation).

Standard values from building codes: - Interior surfaces (still air): \(R_{si} \approx 0.12\) m\(^2\)\(\cdot\)K/W \(\rightarrow\) \(h \approx 8\) W/m\(^2\)\(\cdot\)K - Exterior surfaces (with wind): \(R_{se} \approx 0.03-0.06\) m\(^2\)\(\cdot\)K/W \(\rightarrow\) \(h \approx 17-33\) W/m\(^2\)\(\cdot\)K

These resistances connect the solid surface node to the adjacent air node.

8.6.3 Radiation Heat Transfer

Radiation between surfaces is inherently nonlinear (\(\dot{Q} \propto T^4\)), but for small temperature differences it can be linearized:

\[\dot{Q}_{rad} = h_r A (T_1 - T_2)\]

where the linearized radiation coefficient is:

\[h_r = 4\varepsilon\sigma T_m^3\]

with \(T_m\) being the mean absolute temperature and \(\sigma = 5.67 \times 10^{-8}\) W/m\(^2\)\(\cdot\)K\(^4\).

At typical building temperatures (~300 K) with high emissivity surfaces (\(\varepsilon \approx 0.9\)):

\[h_r \approx 4 \times 0.9 \times 5.67 \times 10^{-8} \times 300^3 \approx 5.5 \text{ W/m}^2\cdot\text{K}\]

The corresponding radiation resistance is:

\[R_{rad} = \frac{1}{h_r A}\]

In networks, radiation often appears in parallel with convection at surfaces, since both transfer heat between the same nodes.

8.6.4 Lumped Building Mass (Thermal Capacitance)

The thermal capacitance at a node represents the heat storage capacity of the associated mass:

\[C = mc = \rho V c\]

where \(m\) is mass, \(V\) is volume, \(\rho\) is density, and \(c\) is specific heat.

Where to place capacitances matters for accuracy:

• At the midpoint of a material layer: Good for symmetric boundary conditions
• At surfaces: Good for capturing surface temperature dynamics
• Split between surface and interior: The common “3R2C” wall model

8.6.4.1 Common Wall Models

2R1C (simplest): One capacitance at the wall midpoint, resistances on either side

T_out ---[R/2]---o---[R/2]--- T_in
                 |
                === C
                 |
                GND

3R2C (more accurate): Splits the capacitance between interior and exterior portions

T_out ---[R1]---o---[R2]---o---[R3]--- T_in
                |          |
               === C1     === C2
                |          |
               GND        GND

The 3R2C model captures the fact that the inner portion of a wall responds differently than the outer portion to temperature changes.

8.6.5 Solar Radiation as a Heat Source

Solar radiation absorbed by a surface enters the network as a heat source (current source in the electrical analogy):

\[\dot{Q}_{sol} = \alpha_s \cdot I_{sol} \cdot A\]

where:

• \(\alpha_s\) is the solar absorptance of the surface
• \(I_{sol}\) is the incident solar irradiance (W/m\(^2\))
• \(A\) is the surface area

In the network, this appears as a current injected at the surface node:

      I_sol (time-varying)
           |
           v
           o---[R]---o---[R]--- ...
           |
          === C
           |
          GND

Solar gains are time-varying inputs that drive the system—they’re known disturbances (from weather data) rather than states to be solved for.

8.6.6 HVAC Heat Extraction/Addition

HVAC systems add or remove heat from the building. They appear as controlled heat sources/sinks:

\[\dot{Q}_{HVAC}(t)\]

For air-based systems (forced air heating/cooling), the heat is typically injected at the indoor air node.

For radiant systems (radiant floors, ceiling panels), the heat is injected at the appropriate surface or mass node.

In control terms, \(\dot{Q}_{HVAC}\) is the control input—it’s what we manipulate to achieve desired temperatures.

8.6.7 Infiltration and Ventilation

Air exchange with outdoors brings in outdoor air at outdoor temperature. The heat flow is:

\[\dot{Q}_{inf} = \dot{m} c_p (T_{out} - T_{in})\]

where \(\dot{m}\) is the mass flow rate of infiltrating air (kg/s) and \(c_p \approx 1005\) J/kg\(\cdot\)K.

This can be represented as a thermal resistance connecting indoor and outdoor air:

\[R_{inf} = \frac{1}{\dot{m} c_p}\]

Note that \(\dot{m}\) depends on wind speed, stack effect, and mechanical ventilation—it may be time-varying.

In the network:

T_out ----[R_inf]---- T_in

Infiltration acts in parallel with the envelope heat loss path.

TipLearn-by-Doing Activity: Drawing a Complete Building Network

Draw the thermal resistance-capacitance network for a simple single-zone building with:

• One exterior wall (modeled as 3R2C)
• One window (resistance only, negligible thermal mass)
• Infiltration (at 0.5 air changes per hour for a 100 m\(^3\) room)
• Internal gains from occupants and equipment: 500 W
• HVAC heating/cooling system

Tasks:

1. Sketch the network showing all R and C elements, heat source inputs, and node temperatures.

2. Label each element with its physical meaning (e.g., “exterior surface resistance,” “wall thermal mass”).

3. Identify which quantities are:

   • States (temperatures at capacitance nodes)
   • Inputs (outdoor temperature, solar radiation, HVAC power)
   • Parameters (resistances, capacitances)

4. Write the KCL equation for the indoor air node.

Hint: Your network should have at least 3 capacitance nodes (2 in the wall, 1 for indoor air) and at least 5 resistance elements.

8.7 Solving the 1R1C Network

8.7.1 The Simplest Building Model

The 1R1C model reduces a building to its bare essentials: a single thermal mass (the building’s interior and structure) connected to the outdoor environment through a single thermal resistance (the envelope).

T_out --------[R]-------- T_in
                           |
                          === C
                           |
                          GND

We’ve already seen this model in action in our motivating example. Now let’s develop the mathematical framework to understand why it behaves as it does.

Despite its simplicity, the 1R1C model captures the essential first-order dynamics of buildings and is widely used for:

• Quick energy estimates
• Understanding fundamental building behavior
• Initial control system design
• Benchmarking more complex models

8.7.2 Deriving the Governing Equation

Applying Kirchhoff’s Current Law (energy conservation) at the indoor node:

\[\underbrace{C\frac{dT_{in}}{dt}}_{\text{rate of energy storage}} = \underbrace{\frac{T_{out} - T_{in}}{R}}_{\text{heat flow through envelope}} + \underbrace{\dot{Q}_{gains}}_{\text{internal + HVAC gains}}\]

Rearranging into standard form:

\[\frac{dT_{in}}{dt} = -\frac{1}{RC}T_{in} + \frac{1}{RC}T_{out} + \frac{1}{C}\dot{Q}_{gains}\]

This is a first-order linear ordinary differential equation (ODE). The coefficient \(-1/(RC)\) determines the system’s dynamics.

8.7.3 The Time Constant

The time constant emerges naturally from the equation:

\[\tau = RC\]

It has units of time (check: [K/W] × [J/K] = [J/W] = [s]) and completely characterizes the building’s dynamic response.

Physical meaning of \(\tau\):

• \(\tau\) is the time required for the system to complete 63.2% of its response to a step change (since \(1 - e^{-1} \approx 0.632\))
• After \(2\tau\): 86.5% complete
• After \(3\tau\): 95.0% complete
• After \(4\tau\): 98.2% complete (essentially at steady state)

What determines \(\tau\)?

\[\tau = RC = \frac{1}{UA} \times mc = \frac{mc}{UA}\]

• Large thermal mass (\(m\) or \(c\)) \(\rightarrow\) large \(\tau\) \(\rightarrow\) slow response
• High envelope conductance (\(UA\)) \(\rightarrow\) small \(R\) \(\rightarrow\) small \(\tau\) \(\rightarrow\) fast response

8.7.4 Analytical Solution

For constant outdoor temperature \(T_{out}\) and constant heat gains \(\dot{Q}\), the ODE has an analytical solution:

\[\boxed{T_{in}(t) = T_{in,\infty} + (T_{in,0} - T_{in,\infty})e^{-t/\tau}}\]

where:

• \(T_{in,0}\) is the initial indoor temperature at \(t=0\)
• \(T_{in,\infty}\) is the steady-state temperature as \(t \to \infty\)

The steady-state temperature is found by setting \(dT_{in}/dt = 0\):

\[T_{in,\infty} = T_{out} + R \cdot \dot{Q}_{gains}\]

This makes physical sense: at steady state, the indoor temperature exceeds outdoor by an amount proportional to the heat gains and the thermal resistance.

Connection to the Fourier heat equation: Recall from our earlier analysis (Section 8.3) that the general solution to the transient heat equation has the form:

\[T(x,t) = T_{ss}(x) + \sum_{n=1}^{\infty} A_n \phi_n(x) e^{-t/\tau_n}\]

This was a steady-state solution plus an infinite sum of exponentially decaying modes, each with its own spatial pattern \(\phi_n(x)\) and time constant \(\tau_n\).

The 1R1C solution above is exactly this structure, but with only a single mode (\(n=1\))! By lumping the building into one thermal mass, we’ve collapsed the infinite series down to a single exponential term:

• \(T_{in,\infty}\) plays the role of \(T_{ss}\) (the steady-state solution)
• \((T_{in,0} - T_{in,\infty})\) is the coefficient \(A_1\) (determined by initial conditions)
• \(e^{-t/\tau}\) is the single exponential decay term
• The spatial variation \(\phi_1(x)\) has disappeared—we have only one temperature, not a temperature distribution

We’ve lost the spatial detail but kept the essential dynamic behavior: exponential decay toward steady state with a characteristic time constant. This is the power of lumped parameter models—they capture the dominant physics while remaining simple enough to analyze and use for control.

Special cases:

1. Free cooling (no heating, \(\dot{Q} = 0\)): \(T_{in,\infty} = T_{out}\) \[T_{in}(t) = T_{out} + (T_{in,0} - T_{out})e^{-t/\tau}\]

2. Constant heating to maintain setpoint: Solve for required \(\dot{Q}\): \[\dot{Q}_{required} = \frac{T_{setpoint} - T_{out}}{R} = UA \cdot (T_{setpoint} - T_{out})\]

8.7.5 Physical Interpretation

The time constant tells us how a building will behave:

Building Type	Typical \(\tau\)	Characteristics
Lightweight office (steel frame, curtain wall)	2-5 hours	Responds quickly to heating/cooling; loses heat fast when HVAC stops
Standard residential	5-15 hours	Moderate response; can coast through short HVAC outages
Heavy masonry/concrete	15-50 hours	Very slow response; excellent thermal coasting; slow to warm up
Historic stone building	50-100+ hours	Extremely stable; seasonal rather than daily dynamics dominate

Example: Comparing two buildings

Consider two buildings with the same envelope resistance (\(R = 0.01\) K/W, i.e., \(UA = 100\) W/K) but different thermal mass:

• Lightweight: \(C = 500\) kJ/K \(\rightarrow\) \(\tau = 500,000/100 = 5,000\) s \(\approx\) 1.4 hours
• Heavyweight: \(C = 5,000\) kJ/K \(\rightarrow\) \(\tau = 5,000,000/100 = 50,000\) s \(\approx\) 14 hours

If both start at 20\(^\circ\)C and outdoor temperature drops to 0\(^\circ\)C with no heating:

• After 4 hours: Lightweight at 11.4\(^\circ\)C, Heavyweight at 16.7\(^\circ\)C
• After 8 hours: Lightweight at 6.5\(^\circ\)C, Heavyweight at 14.0\(^\circ\)C

The heavyweight building “coasts” much longer—critical for resilience during power outages or for taking advantage of thermal mass for load shifting.

TipLearn-by-Doing Activity: 1R1C Model Analysis

Return to our motivating example: the office room with \(UA = 7.08\) W/K and \(C = 52\) kJ/K (air only).

Tasks:

1. Verify that \(\tau \approx 2\) hours using \(\tau = C/UA\).

2. If heating fails at midnight when \(T_{in} = 20^\circ\text{C}\) and \(T_{out} = 0^\circ\text{C}\), write the equation for \(T_{in}(t)\) and calculate:

   • Temperature at 2 AM (after 2 hours)
   • Temperature at 6 AM (after 6 hours)

3. How long until \(T_{in}\) drops to 15\(^\circ\)C? Solve \(15 = 0 + 20 \cdot e^{-t/\tau}\) for \(t\).

4. What heating power is required to maintain 20\(^\circ\)C at steady state when \(T_{out} = 0^\circ\text{C}\)?

5. If we add the concrete floor (\(C_{total} = 4,120\) kJ/K), recalculate \(\tau\) and the time to reach 15\(^\circ\)C.

Answers:

1. \(\tau = 52,000/7.08 = 7,345\) s \(\approx\) 2.04 hours \(\checkmark\)

2. \(T_{in}(t) = 0 + 20 \cdot e^{-t/7345}\)

   • At 2 AM (\(t = 7200\) s): \(T_{in} = 20 \cdot e^{-0.98} = 7.5^\circ\text{C}\)
   • At 6 AM (\(t = 21600\) s): \(T_{in} = 20 \cdot e^{-2.94} = 1.1^\circ\text{C}\)

3. \(15 = 20 \cdot e^{-t/\tau}\) \(\rightarrow\) \(t = -\tau \ln(0.75) = 0.288\tau = 35\) minutes

4. \(\dot{Q} = UA \cdot \Delta T = 7.08 \times 20 = 142\) W

5. With concrete: \(\tau = 4,120,000/7.08 = 582,000\) s \(\approx\) 162 hours (!), time to 15\(^\circ\)C \(\approx\) 47 hours

The massive increase in time constant with concrete dramatically changes the building’s behavior!

8.8 Solving the 2R2C Network

8.8.1 Why Two Capacitances?

The 1R1C model lumps all thermal mass into a single capacitance. But real buildings have distinct thermal masses that respond at different rates:

• Indoor air: Low thermal mass, responds in minutes to hours
• Furnishings and interior surfaces: Medium mass, responds in hours
• Building structure (walls, floors): High thermal mass, responds in hours to days

A 2R2C model captures this separation of time scales with two capacitances, enabling more realistic predictions of: - Short-term temperature swings (air heats up fast, but walls lag) - The “warm up” period after turning on heating - How quickly a room “feels” comfortable versus how quickly it reaches thermal equilibrium

8.8.2 Network Configuration

The most common 2R2C configuration for buildings separates the indoor air from the building mass:

T_out ----[R_1]---- T_wall ----[R_2]---- T_air
                       |                   |
                      === C_wall          === C_air
                       |                   |
                      GND                 GND

Where:

• \(T_{air}\) = indoor air temperature (what occupants feel immediately)
• \(T_{wall}\) = average temperature of interior surfaces and structure
• \(C_{air}\) = thermal capacitance of air + furnishings (small, ~50-200 kJ/K for a room)
• \(C_{wall}\) = thermal capacitance of building mass (large, ~500-5000 kJ/K)
• \(R_1\) = resistance from outdoor to wall mass (exterior envelope)
• \(R_2\) = resistance from wall mass to indoor air (interior surface film)

Heat gains (HVAC, occupants, solar) typically enter at the air node or wall node depending on whether they’re convective or radiative.

8.8.3 Deriving the System of ODEs

Applying KCL at each capacitance node:

At the wall node (\(T_w\)): \[C_w \frac{dT_w}{dt} = \frac{T_{out} - T_w}{R_1} + \frac{T_a - T_w}{R_2} + \dot{Q}_{w}\]

At the air node (\(T_a\)): \[C_a \frac{dT_a}{dt} = \frac{T_w - T_a}{R_2} + \dot{Q}_{a}\]

Where \(\dot{Q}_w\) represents heat gains absorbed by the walls (e.g., solar on interior surfaces) and \(\dot{Q}_a\) represents gains to the air (e.g., HVAC, occupants).

These are two coupled first-order ODEs—the temperature of each node depends on the temperature of the other.

8.8.4 Two Time Constants

Unlike the 1R1C model with its single time constant, the 2R2C system has two time constants: \(\tau_{fast}\) and \(\tau_{slow}\).

Intuitively:

• \(\tau_{fast}\): Dominated by the smaller capacitance (air). Governs the rapid initial response when heating starts.
• \(\tau_{slow}\): Dominated by the larger capacitance (structure). Governs the slow approach to final steady state.

Typically \(\tau_{slow} / \tau_{fast} \approx 10-100\), meaning the system has two distinct response phases:

1. Fast phase (minutes to ~1 hour): Air temperature changes rapidly; walls barely respond
2. Slow phase (hours to days): Air and wall temperatures gradually equilibrate

This two-phase behavior explains why:

• A room can feel warm (air heated) even though walls are still cold
• Radiant heating (heats walls directly) feels different than forced-air heating (heats air)
• Buildings need “pre-conditioning” time before occupancy

8.8.5 Solving the System

To find the time constants, we write the system in matrix form. Let \(\mathbf{x} = [T_w, T_a]^T\) be the state vector. The ODEs become:

\[\frac{d\mathbf{x}}{dt} = \mathbf{A}\mathbf{x} + \mathbf{B}\mathbf{u}\]

where \(\mathbf{u}\) contains the inputs (\(T_{out}\), \(\dot{Q}_w\), \(\dot{Q}_a\)) and:

\[\mathbf{A} = \begin{bmatrix} -\frac{1}{C_w}\left(\frac{1}{R_1} + \frac{1}{R_2}\right) & \frac{1}{C_w R_2} \\ \frac{1}{C_a R_2} & -\frac{1}{C_a R_2} \end{bmatrix}\]

The eigenvalues of \(\mathbf{A}\) determine the time constants:

\[\lambda_1, \lambda_2 = \text{eigenvalues of } \mathbf{A}\] \[\tau_1 = -\frac{1}{\lambda_1}, \quad \tau_2 = -\frac{1}{\lambda_2}\]

For thermal systems, both eigenvalues are negative (stable system) and real (no oscillations).

The general solution involves the matrix exponential:

\[\mathbf{x}(t) = e^{\mathbf{A}t}\mathbf{x}(0) + \int_0^t e^{\mathbf{A}(t-s)}\mathbf{B}\mathbf{u}(s)\,ds\]

For constant inputs, this simplifies to exponential decay toward steady state—but now with two exponential terms, each with its own time constant.

8.8.6 Physical Interpretation of Multiple Time Constants

The two-time-constant behavior has important practical implications:

When heating starts:

1. Air temperature rises quickly (fast time constant)
2. Walls remain cold, absorbing heat from the warm air
3. Gradually, walls warm up (slow time constant)
4. Eventually, everything equilibrates

When heating stops:

1. Air temperature drops quickly (fast time constant)
2. But warm walls release stored heat, slowing the air temperature drop
3. Room stays warmer longer than a 1R1C model would predict

This is why occupants might feel cold shortly after heating starts (walls radiating cold toward them) even though the air thermostat shows a comfortable temperature.

TipLearn-by-Doing Activity: 2R2C Model Simulation

Consider a room with the following 2R2C parameters:

• \(C_{air} = 50\) kJ/K (air and light furnishings)
• \(C_{wall} = 500\) kJ/K (interior wall mass)
• \(R_1 = 0.05\) K/W (exterior envelope resistance)
• \(R_2 = 0.01\) K/W (interior surface resistance)
• Initial conditions: \(T_{air} = T_{wall} = 15^\circ\text{C}\)
• Outdoor temperature: \(T_{out} = 0^\circ\text{C}\)
• HVAC heating: \(\dot{Q} = 1000\) W applied to air node at \(t=0\)

Tasks:

1. Write the system of ODEs in matrix form: \(\dot{\mathbf{x}} = \mathbf{A}\mathbf{x} + \mathbf{B}\mathbf{u}\)

2. Calculate the eigenvalues of \(\mathbf{A}\) and determine the two time constants.

3. Using Python (or MATLAB), simulate the temperature response for 24 hours. Plot both \(T_{air}(t)\) and \(T_{wall}(t)\).

4. Create a 1R1C model with equivalent total R (\(R_{total} = R_1 + R_2\)) and C (\(C_{total} = C_{air} + C_{wall}\)). Simulate and compare.

5. What’s different between the 1R1C and 2R2C responses? Which is more realistic?

Hints:

• The A matrix should be 2×2
• Use numpy.linalg.eig() for eigenvalues
• Use scipy.integrate.solve_ivp() for simulation
• Time constants: \(\tau = -1/\lambda\) for each eigenvalue

Expected results:

• Fast time constant: ~10-20 minutes
• Slow time constant: ~5-10 hours
• Air temperature overshoots initially then settles as walls warm up
• 1R1C model misses this overshoot behavior

8.9 State-Space Representation

8.9.1 Why State-Space?

State-space representation is the standard language for describing dynamic systems in modern control theory. Converting our thermal network models to state-space form unlocks powerful tools:

• Stability analysis: Are all eigenvalues of A negative? (For thermal systems, yes)
• Controllability: Can we drive temperatures to any desired state using available inputs?
• Observability: Can we estimate all temperatures from available sensor measurements?
• State feedback control: Design controllers that use full state information
• State estimation: Kalman filters to estimate unmeasured temperatures
• Model Predictive Control (MPC): Optimize control actions over a prediction horizon

For your final projects involving autonomous building control, state-space models will be essential.

8.9.2 The Standard Form

A linear time-invariant (LTI) system in state-space form consists of two equations:

\[\boxed{\dot{\mathbf{x}} = \mathbf{A}\mathbf{x} + \mathbf{B}\mathbf{u}}\] (State equation)

\[\boxed{\mathbf{y} = \mathbf{C}\mathbf{x} + \mathbf{D}\mathbf{u}}\] (Output equation)

Where:

Symbol	Name	Description	Dimensions
\(\mathbf{x}\)	State vector	Temperatures at capacitance nodes	\(n \times 1\)
\(\mathbf{u}\)	Input vector	External inputs (weather, HVAC, gains)	\(m \times 1\)
\(\mathbf{y}\)	Output vector	Measured quantities	\(p \times 1\)
\(\mathbf{A}\)	System matrix	Describes internal dynamics	\(n \times n\)
\(\mathbf{B}\)	Input matrix	How inputs affect states	\(n \times m\)
\(\mathbf{C}\)	Output matrix	What we measure	\(p \times n\)
\(\mathbf{D}\)	Feedthrough matrix	Direct input-to-output path	\(p \times m\)

For thermal systems, \(\mathbf{D} = \mathbf{0}\) (no instantaneous effect of inputs on outputs).

8.9.3 Converting Thermal Networks to State-Space

Systematic procedure:

1. Identify state variables: One state per capacitance node (the temperatures)
2. Write KCL equations: Energy balance at each capacitance node
3. Identify inputs: Outdoor temperature, solar radiation, HVAC power, internal gains
4. Rearrange into matrix form: Isolate \(d\mathbf{x}/dt\) on the left
5. Extract A, B, C, D matrices: Read off coefficients

Key insight: The structure of \(\mathbf{A}\) comes directly from the network topology. Each entry \(A_{ij}\) represents how temperature at node \(j\) affects the rate of change of temperature at node \(i\).

8.9.4 Worked Example: 1R1C in State-Space

Starting from our 1R1C equation:

\[C\frac{dT_{in}}{dt} = \frac{T_{out} - T_{in}}{R} + \dot{Q}\]

Divide by \(C\) and rearrange:

\[\frac{dT_{in}}{dt} = -\frac{1}{RC}T_{in} + \frac{1}{RC}T_{out} + \frac{1}{C}\dot{Q}\]

Define state and inputs:

• State: \(x = T_{in}\) (scalar, \(n=1\))
• Inputs: \(\mathbf{u} = \begin{bmatrix} T_{out} \\ \dot{Q} \end{bmatrix}\) (\(m=2\))

State-space matrices:

\[A = -\frac{1}{RC} = -\frac{1}{\tau}\]

\[B = \begin{bmatrix} \frac{1}{RC} & \frac{1}{C} \end{bmatrix} = \begin{bmatrix} \frac{1}{\tau} & \frac{1}{C} \end{bmatrix}\]

If we measure indoor temperature directly: \(C = 1\), \(D = \begin{bmatrix} 0 & 0 \end{bmatrix}\)

The single eigenvalue of \(A\) is \(\lambda = -1/\tau\), confirming our time constant.

8.9.5 Worked Example: 2R2C in State-Space

For the 2R2C network with states \(\mathbf{x} = \begin{bmatrix} T_w \\ T_a \end{bmatrix}\) and inputs \(\mathbf{u} = \begin{bmatrix} T_{out} \\ \dot{Q}_a \end{bmatrix}\):

Starting from the nodal equations:

\[C_w \frac{dT_w}{dt} = \frac{T_{out} - T_w}{R_1} + \frac{T_a - T_w}{R_2}\]

\[C_a \frac{dT_a}{dt} = \frac{T_w - T_a}{R_2} + \dot{Q}_a\]

Dividing by capacitances and rearranging:

\[\frac{dT_w}{dt} = -\frac{1}{C_w}\left(\frac{1}{R_1} + \frac{1}{R_2}\right)T_w + \frac{1}{C_w R_2}T_a + \frac{1}{C_w R_1}T_{out}\]

\[\frac{dT_a}{dt} = \frac{1}{C_a R_2}T_w - \frac{1}{C_a R_2}T_a + \frac{1}{C_a}\dot{Q}_a\]

State-space matrices:

\[\mathbf{A} = \begin{bmatrix} -\frac{1}{C_w}\left(\frac{1}{R_1} + \frac{1}{R_2}\right) & \frac{1}{C_w R_2} \\ \frac{1}{C_a R_2} & -\frac{1}{C_a R_2} \end{bmatrix}\]

\[\mathbf{B} = \begin{bmatrix} \frac{1}{C_w R_1} & 0 \\ 0 & \frac{1}{C_a} \end{bmatrix}\]

The eigenvalues of \(\mathbf{A}\) give the two time constants: \(\tau_i = -1/\lambda_i\).

8.9.6 Properties of the System Matrix A

For thermal systems, the system matrix \(\mathbf{A}\) has special properties:

1. Stability (all eigenvalues have negative real parts)

Physical systems that dissipate energy are inherently stable. Heat always flows from hot to cold, so temperature differences decay over time. Mathematically, this means all eigenvalues of \(\mathbf{A}\) satisfy \(\text{Re}(\lambda_i) < 0\).

2. Real eigenvalues (no oscillations)

Pure thermal systems don’t oscillate—temperatures decay monotonically toward equilibrium. This is because \(\mathbf{A}\) for thermal networks is a special type of matrix (related to graph Laplacians) that guarantees real eigenvalues.

Note: Coupled thermo-fluid systems or systems with feedback control CAN oscillate, but the thermal subsystem alone cannot.

3. Eigenvalues give time constants

Each eigenvalue \(\lambda_i\) corresponds to a time constant:

\[\tau_i = -\frac{1}{\lambda_i}\]

The eigenvalue with smallest magnitude (closest to zero) gives the slowest time constant, which dominates the long-term response.

4. Structure reflects network topology

• Diagonal entries \(A_{ii}\) are negative (self-damping)
• Off-diagonal entries \(A_{ij}\) are positive if nodes \(i\) and \(j\) are connected
• Row sums relate to heat loss to external nodes (boundary conditions)

TipLearn-by-Doing Activity: State-Space Formulation

Consider the following 3-node thermal network representing a room with a wall:

T_out ---[R_1]--- T_wall ---[R_2]--- T_air ---[R_3]--- T_adj
                    |                  |
                   === C_w            === C_a
                    |                  |
                   GND                GND

Where:

• \(T_{out}\) = outdoor temperature (input)
• \(T_{adj}\) = adjacent room temperature (input, assume constant at 20\(^\circ\)C)
• \(T_{wall}\) = wall mass temperature (state)
• \(T_{air}\) = room air temperature (state)
• \(\dot{Q}_{hvac}\) = HVAC heat input to air node (input)

Parameters:

• \(R_1 = 0.05\) K/W, \(R_2 = 0.02\) K/W, \(R_3 = 0.1\) K/W
• \(C_w = 800\) kJ/K, \(C_a = 100\) kJ/K

Tasks:

1. Write the KCL equation at each capacitance node.

2. Define \(\mathbf{x} = \begin{bmatrix} T_w \\ T_a \end{bmatrix}\) and \(\mathbf{u} = \begin{bmatrix} T_{out} \\ T_{adj} \\ \dot{Q}_{hvac} \end{bmatrix}\).

3. Express the system in the form \(\dot{\mathbf{x}} = \mathbf{A}\mathbf{x} + \mathbf{B}\mathbf{u}\).

4. Calculate the numerical values of \(\mathbf{A}\) and \(\mathbf{B}\).

5. Find the eigenvalues of \(\mathbf{A}\) and the corresponding time constants.

6. If we measure only \(T_{air}\), what is \(\mathbf{C}\)?

Python verification:

import numpy as np

# Parameters (convert kJ to J)
R1, R2, R3 = 0.05, 0.02, 0.1  # K/W
Cw, Ca = 800e3, 100e3  # J/K

# A matrix
A = np.array([
    [-(1/Cw)*(1/R1 + 1/R2), (1/Cw)*(1/R2)],
    [(1/Ca)*(1/R2), -(1/Ca)*(1/R2 + 1/R3)]
])

# Eigenvalues and time constants
eigenvalues = np.linalg.eigvals(A)
time_constants = -1/eigenvalues
print(f"Eigenvalues: {eigenvalues}")
print(f"Time constants: {time_constants/3600} hours")

Eigenvalues: [-3.24397649e-05 -6.55060235e-04]
Time constants: [8.56287889 0.42404921] hours

Expected results:

• Two time constants: one fast, air-dominated, and one slow, wall-dominated.

8.10 Putting It All Together

8.10.1 From Physics to Control

Let’s step back and appreciate the journey we’ve taken:

Fourier Heat Equation (PDE)
    ↓ Spatial discretization
Thermal Network (System of ODEs)
    ↓ Matrix formulation
State-Space Model (ẋ = Ax + Bu)
    ↓ Ready for...
Control Design, Simulation, Identification

1. Physics \(\rightarrow\) PDE: We started with the fundamental law of heat conduction (Fourier’s equation), which describes how temperature varies continuously in space and time.

2. PDE \(\rightarrow\) Network: By discretizing space—either uniformly or at physical boundaries—we converted the infinite-dimensional PDE into a finite-dimensional system of ODEs. The thermal network representation makes the physics intuitive: resistances impede heat flow, capacitances store heat.

3. Network \(\rightarrow\) State-Space: By writing the nodal equations in matrix form, we arrived at the standard state-space representation used throughout control theory. This unlocks decades of tools for analysis and design.

This progression—from first-principles physics to control-ready models—is fundamental not just to buildings, but to most engineered systems. The skills you’ve learned here apply to thermal management in electronics, chemical process control, vehicle cabin climate, and beyond.

8.10.2 Limitations and Extensions

Our models make simplifying assumptions. It’s important to understand when they apply and when more sophisticated approaches are needed:

Linearity: We assumed constant R and C, independent of temperature. This is excellent for typical building temperature ranges ($\(20\)^$C from nominal) but breaks down for:

• High-temperature processes
• Phase-change materials
• Systems with strong radiation coupling

Single-zone: We modeled single rooms. Real buildings have multiple zones with:

• Inter-zone heat transfer
• Shared HVAC systems
• Varying occupancy patterns

No humidity: We focused on sensible heat (temperature). Latent heat (humidity) matters for:

• Comfort (humid air feels warmer)
• Cooling loads (dehumidification is energy-intensive)
• Condensation risk

Known parameters: We assumed R and C are known. In practice:

• Material properties vary
• Construction details are uncertain
• Parameter identification from data is often needed (a topic for later in this course)

Neglected dynamics: We ignored:

• HVAC system dynamics (heating/cooling isn’t instantaneous)
• Sensor and actuator dynamics
• Air mixing and stratification

8.10.3 Preview of What’s Next

In the following lectures, we will use these thermal models to:

• Understand occupant thermal comfort (what temperatures should we target?)
• Design and analyze HVAC control strategies
• Apply model predictive control to optimize energy use while maintaining comfort

8.11 Additional Resources

8.11.1 Primary Reference

• Chapter 8 of Reddy (as noted in syllabus)

8.11.2 Supplementary Reading

• Lecture 5 notes (prerequisite steady-state material)
• Any introductory control systems textbook for state-space fundamentals

8.11.3 Online Resources

8.11.3.1 Thermal Modeling

• Engineering ToolBox - Thermal Diffusivity - Material property data
• Thermtest - Heat Penetration Calculator - Online calculator for penetration depth

8.11.3.2 State-Space and Control

• MIT OpenCourseWare - Introduction to State-Space - Free course materials
• Python Control Systems Library - Tools for control system analysis in Python

8.11.3.3 Building Simulation

• EnergyPlus - Detailed building energy simulation (uses thermal networks internally)
• Modelica Buildings Library - Component-based thermal modeling