Skip to content

Instantly share code, notes, and snippets.

@ed-davies
Last active Jan 3, 2016
Embed
What would you like to do?
Quick-and-dirty simulation of the effect of thermal mass on a cartoon house's heating requirements.
#!/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()
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment