Quick-and-dirty simulation of the effect of thermal mass on a cartoon house's heating requirements.

thermal-mass.py

#!/usr/bin/env python3

import math

# Constants to reference times of day in seconds conveniently.
MINUTES = 60
HOURS = 60 * MINUTES
DAYS = 24 * HOURS

# House dimensions (in metres) and similar parameters.
class house:
    eastWest = 7
    northSouth = 6
    height = 4.8
    u_wall = 0.1                # Heat loss through the walls, roof and floor.
    u_window = 1.0              # Heat loss through the windows.
    glazing = 0.5               # Proportion of south face glazed.
    thermalMass = 20            # Thermal mass expressed as an equivalent amount
                                # in millimetres of 2 MJ/(m³·K) material (like 
                                # concrete) on the floor.

# Simulation parameters.
class sim:
    timeStep = 5 * MINUTES      # Time increment over which we assume stable
                                # flows depending on the start conditions.
    trainingPeriod = 50 * DAYS  # Period over which we run the simulation to
                                # reach initial near-stable conditions.
    runLength = 1 * DAYS        # Period over which we collect data.

# Utility functions.
def pairs(it):
    """ Generate consecutive pairs in an iterable. """
    it = iter(it)
    p = next(it)
    for n in it:
        yield p, n
        p = n
        
def interpolate(tod, values):
    """ Find the value of a variable at a specified time of day by
        linearly interpolating between appropriate entries in a list
        of (time of day, value) pairs.
    """
    for (start, startValue), (end, endValue) in pairs(values):
        if (start <= tod) and (tod <= end):
            fraction = (tod - start) / (end - start)
            return startValue + (endValue - startValue) * fraction
            
def timeString(t):
    """ Format a time in seconds in a reasonably readable form. """
    if t >= DAYS:
        d = math.floor(t / DAYS)
        return "%dd%s" % (d, timeString(t - (d * DAYS)))
    h = math.floor(t / HOURS)
    t -= h * HOURS
    m = math.floor(t / MINUTES)
    t -= m * MINUTES
    if t != 0:
        return "%02d:%02d:%02f" % (h, m, t)
    else:
        return "%02d:%02d" % (h, m)

# Environmental conditions as a function of time of day (in seconds).
def outsideTemperature(tod):
    night = 2.0                 # Temperature over night.
    afternoon = 5.0             # Temperature during the afternoon.
    return interpolate(tod, (
        (0 * HOURS, night),
        (6 * HOURS, night),
        (12 * HOURS, afternoon),
        (18 * HOURS, afternoon),
        (24 * HOURS, night)))
        
def insolation(tod):
    """ Sunlight on a vertical south facing surface in W/m². """
    return interpolate(tod, (
        (0 * HOURS, 0),
        (10 * HOURS, 0),
        (12 * HOURS, 300),
        (14 * HOURS, 0),
        (24 * HOURS, 0)))
        
def incidentalGain(tod):
    """ Heating from bodies, lights, appliances, etc. """
    return interpolate(tod, (
        (0 * HOURS, 100),
        (7 * HOURS, 100),
        (8 * HOURS, 200),
        (22 * HOURS, 200),
        (23 * HOURS, 100),
        (24 * HOURS, 100)))
        
def tempRequired(tod):
    """ Temperature required in the house during the 24 hour period. """
    return interpolate(tod, (
        (0 * HOURS, 15),
        (7 * HOURS, 15),
        (8 * HOURS, 19),
        (22 * HOURS, 19), 
        (23 * HOURS, 15),
        (24 * HOURS, 15)))
        
# Simulation.
def timeOfSim():
    """ Yield the times of simulated day, together with whether or not
        we're past the end of the training period. """
    simLength = sim.trainingPeriod + sim.runLength
    tos = 0
    while tos < simLength:
        yield tos % DAYS, tos >= sim.trainingPeriod
        tos += sim.timeStep

def runSimulation():
    # Calculate derived house parameters.
    heatCapacity = (house.eastWest * house.northSouth * house.thermalMass /
                    1000 * 2e6)
    windowArea = house.eastWest * house.height * house.glazing
    wallArea = (((house.eastWest + house.northSouth) * 2 * house.height) +
                (house.eastWest * house.northSouth * 2) -
                windowArea)
    heatLossCoefficient = wallArea * house.u_wall + windowArea * house.u_window
    
    # Run simulation.
    temp = 19
    
    avg = 0
    count = 0
    min_temp = 10000
    max_temp = -10000
    heating = 0
    
    for tod, collectingData in timeOfSim():

        heatLoss = (heatLossCoefficient * (temp - outsideTemperature(tod)) -
                    windowArea * insolation(tod) -
                    incidentalGain(tod)) * sim.timeStep
        
        temp -= heatLoss / heatCapacity
        
        tReq = tempRequired(tod)
        
        if temp < tReq:
            heat_required = (tReq - temp) * heatCapacity
            temp = tReq
        else:
            heat_required = 0
        
        if collectingData:
            if (tod % HOURS) == 0:
                print(timeString(tod), temp)
            avg += temp
            count += 1
            min_temp = min(min_temp, temp)
            max_temp = max(max_temp, temp)
            heating += heat_required

    print(avg/count, min_temp, max_temp, heating / 3.6e6)
    
for house.thermalMass in (50, 100, 200, 400):
    print(house.thermalMass)
    runSimulation()
    print()

tm.txt

50
00:00 21.59208674351369
01:00 21.080520577513244
02:00 20.5845185326866
03:00 20.103607079148446
04:00 19.63732709390159
05:00 19.18523342251524
06:00 18.746894454139014
07:00 18.330170560594304
08:00 19.0
09:00 19.0
10:00 19.0
11:00 19.888898209653387
12:00 22.88675390706115
13:00 25.611968301553468
14:00 26.1245577719312
15:00 25.650879383570793
16:00 25.19161240144448
17:00 24.746318366389364
18:00 24.314572159123077
19:00 23.887682743237818
20:00 23.458568977765204
21:00 23.02729853732645
22:00 22.59393703757938
23:00 22.112554480375714
21.89063685814395 18.0625807114554 26.200127053911718 1.9095102454363961

100
00:00 20.738070204709228
01:00 20.493559552283138
02:00 20.252794787310144
03:00 20.01571852303563
04:00 19.782274251866383
05:00 19.552406331901853
06:00 19.326059973671835
07:00 19.107340162698133
08:00 19.0
09:00 19.0
10:00 19.0
11:00 19.446246437448693
12:00 20.96118707353239
13:00 22.364313468222882
14:00 22.673548563618088
15:00 22.4879016909144
16:00 22.30509891634822
17:00 22.125096668519205
18:00 21.947852043538038
19:00 21.76916385917938
20:00 21.585553198491777
21:00 21.397095473495668
22:00 21.20386494090419
23:00 20.982829521928856
20.728576760822968 18.936082682870374 22.69787120903317 0.905490321351111

200
00:00 19.997472618116568
01:00 19.880477780307782
02:00 19.764382289839237
03:00 19.649179233366254
04:00 19.53486175068753
05:00 19.421423034336584
06:00 19.30885632917639
07:00 19.19923930087215
08:00 19.10588826754755
09:00 19.026871505459223
10:00 19.0
11:00 19.223574387862634
12:00 19.985090009550078
13:00 20.69698248854616
14:00 20.865307801543235
15:00 20.786056029117997
16:00 20.707413470501457
17:00 20.629375442625296
18:00 20.551937298420267
19:00 20.473010057664446
20:00 20.39084600173656
21:00 20.305470012256773
22:00 20.21690677957853
23:00 20.11360097717547
19.95112323276314 19.0 20.874823577976606 0.23377047310749788

400
00:00 19.709061173762535
01:00 19.65157078624868
02:00 19.594301755772186
03:00 19.537253230035123
04:00 19.480424360021203
05:00 19.423814299983125
06:00 19.36742220742998
07:00 19.31229065577444
08:00 19.26509328149519
09:00 19.224901426051677
10:00 19.1867894882065
11:00 19.296827231944803
12:00 19.67788964917409
13:00 20.035734247742763
14:00 20.122677836442513
15:00 20.08584124626243
16:00 20.049146489157767
17:00 20.012593019024212
18:00 19.97618029186011
19:00 19.93886435309872
20:00 19.899766927423464
21:00 19.858894874149144
22:00 19.816255026179906
23:00 19.766057453112673
19.678738685468723 19.185582127299575 20.12695549305617 0.0