/*
This formulation evaluates heat flow through a single wall comprising a set of
material layers; that it ignores any internal wall structure, studding, or edge
effects at corners, joins with the slab, windows, etc. Because of this
uniform-in-area assumption, any lateral heat flows can be ignored and so the
computation can be treated as in one distance dimension (though the wall) and
one time dimension (∆x,∆t).
By asssuming that each material layer is a fixed multiple a fixed thickness ∆x
and by using a fixed time step ∆t, the solution can be approximated on a fixed
(∆x, ∆t) mesh using a 1-D implementation of a formulation of the 1-D heat
equation. See the Wikipedia "Heat Equation" article for the maths derivation.
The formulation uses 3 material properties that are specific to each material
in the wall:
rho the density
cp the heat capacity
k the thermal conductivity
It also is expressed in terms of the heat flow, W, at each layer boundary and of
the temperature at the centre of each layer, so T(i) is midpoint between
boundaries i and i+1 and conversely W(i) is calculated midpoint between
temperatures T(i-1) and T(i). Note that a postive W indicates a left->right heat
flow (that is T is decreasing in the direction of increasing i). The basic
equations internal to a given layer are straightforward:
The temperature change is proportional to the net heat flow across the layer:
T(i,n+1) = T(i,n) + (del_t / rho * cp * del_x) (W(i,n) W(i+1,n)) (Eq 1)
The heat flow is proportional to the (decreasing) temperature gradient:
W(i+1,n) = (T(i,n) - T(i+1,n)) * k / del_x (Eq 2)
where i is the ∆x index and n the ∆t one. Also note that the second equation can
be substituted in the first giving the standard form derived in the Wikipedia
article.
These equations are in explicit form -- that is given a starting vector of
T(i,0) and suitable boundary conditions, the subsequent T(i,n) can be computed
at following ∆t intervals without the need to solve implicit equations. This
iterative solution is stable for a small domain of (∆x, ∆t) near (0,0), but the
evaluation can become unstable if ∆x or ∆t are too large. My approach to this
issue is simply to use a small enough mesh size, say (1cm, 1s), since cost of
evaluating this on a modern CPU which has GFLOP performance is still measured in
10s of mSec.
However, some care is required in the treatment of material boundaries.
| | |
| Layer i | Layer i+1 |
| Material A | Material B |
| | |
| T(i) | T(i+1) |
W(i) W(i+1) W(i+2)
* Homogenous Materials (A = B). This gives Eq 1 and 2 above.
* Material boundary (A!=B). The W is evaluated at midpoint between two
different k values, so a composite 2/(1/kA + 1/kB) is used in Eq 2.
* Air-Material boundary. (For internal air gaps see the discussion following.)
For external boundaries, a simple Newtons law of cooling approximation is
used:
W = Q ∆T
where Q is a air surface conduction coefficient, and uniform air
temperature is assumed so
∆T = Tair - Tsurface
The surface temperature of the material is simply a linear extropolation
from the 2 adjacent internal temperatures, so if A is air, the surface
temperature corresponding to W(i+1) is
Tsurface = (3*T(i+1) - T(i+2))/2
and if B is the air then
Tsurface = (3*T(i) - T(i-1))/2
This Newtons law approach + specific heat calculation fails for internal air
gaps. This is because the specific heat per unit volume of air is ~2,000 times
smaller than stone, say, so even if a 5mm grid is used, a 50mm air gap has 200
times smaller a specific heat than an adjacent stone layer, so a transfer of
X watts from stone to air-gap resulting in a -∆T in the stone will create
a 200∆T temperature rise in the air. This ill-conditioning causes the mesh
approximation to become unstable for extremely small perturbations. Not good.
So I have adopted the simple approach of ignoring the air gap as a heat
reservoir (as this isn't material in the overall heat capacity calcs anyway) and
calculating the heat transfer directly from the adjacent surfaces:
| | | |
| Layer i | | Layer i+1 |
| Material A | | Material C |
| | | |
| T(i-1) | | T(i) |
W(i-1) W(i) W(i+i)
^ Air gap "optimised" away
W(i) = Q ∆T
= Q (Tleft_surface - Tright_surface)
Also I need to be careful to maintain any dependency ordering in evaluating
the equations -- for example I can't use W(i+1,n) before it itself is
evaluated, etc..
As previously mentioned a fine mesh is used to ensure numeric stability, but
the ouput is reported as a pair of simple CSV files for T and W on a fixed
(∆x, ∆t) grid for subsequent import into Excel or equivalent for graphing and
presentation.
Lastly, I have could has used any of a dozen or more languages to code this. I
had initially intended to use a simple spreadsheet formulation, but the stable
mesh size proved too small for this to be manageable. So instead I have chosen
to use plain C, though I have deliberately kept the coding simple so that this
could easily be reimplemented in any favoured language.
Terry Ellison -- 18 Oct 2014
This is an original work of the author. It is free code and can be used by
anyone without warranty or constraint.
*/
#include
#include
#include
/* Depths are in integer mm */
#define WALL_DEPTH 470
#define LAYER_DEPTH 5
#define LAYERS_PER_M (1000/LAYER_DEPTH)
#define REPORT_DEPTH 50
#define MAX_LAYERS (2 + (WALL_DEPTH/LAYER_DEPTH))
#define TIME_INT_PER_DAY 86400
#define SECS_PER_DAY 86400.0
#define NUM_DAYS 20
#define Q_FACTOR 10.0
#define TWO_PI 6.28318530717959
#define AIR_INTERNAL -1
#define AIR_EXTERNAL -2
#define REPORT_INTERVAL (TIME_INT_PER_DAY / 48)
#define REPORT_END -1
/* *5deg mean with +/- 5deg daily ripple* */
#define REPORT_DELAY (1 * TIME_INT_PER_DAY)
#define START_COLD 0
#define INTERNAL_TEMP 20.0
#define EXTERNAL_TEMP 5.0
#define EXTERNAL_DELTA 5.0
#define EXTERNAL_START ((78 * TIME_INT_PER_DAY) / 24 )
/* *Heating the wall from cold*
#define REPORT_DELAY 0
#define START_COLD 1
#define INTERNAL_TEMP 20.0
#define EXTERNAL_TEMP 5.0
#define EXTERNAL_DELTA 5.0
#define EXTERNAL_START (NUM_DAYS * TIME_INT_PER_DAY)
*/
#define MAX_INTERNAL_HEAT 10.0
enum layer_type {
AIR = 0,
PLASTER_BOARD,
MINERAL_WOOL,
OSB,
CELLULOSE_FIBRE,
PANELVENT,
STONE};
typedef struct {
enum layer_type type;
double K;
double rho;
double cp;
} properties_type;
properties_type property[] = {
/* material type K rho cp */
{AIR, 0.024, 2.0, 2010.0},
{PLASTER_BOARD, 0.170, 3100.0, 2180.0},
{MINERAL_WOOL, 0.040, 150.0, 1400.0},
{OSB, 0.147, 1700.0, 3400.0},
{CELLULOSE_FIBRE, 0.041, 300.0, 3200.0},
{PANELVENT, 0.100, 1000.0, 3400.0},
{STONE, 1.700, 4800.0, 1600.0}
};
/*
* Define the wall layers. The count parameter is overloaded in the case of
* AIR in that it is the actual thickness, -1 (Internal) or -2 (External).
* For the other layers it is the number of LAYER_DEPTH slices.
*/
typedef struct {
enum layer_type type;
int count;
} layer_type;
layer_type layer[] = {
/* material type number of layers */
{AIR, AIR_INTERNAL},
{PLASTER_BOARD, (10 / LAYER_DEPTH)},
{MINERAL_WOOL, (45 / LAYER_DEPTH)},
{OSB, (10 / LAYER_DEPTH)},
{CELLULOSE_FIBRE, (240 / LAYER_DEPTH)},
{PANELVENT, (10 / LAYER_DEPTH)},
{AIR, (50 / LAYER_DEPTH)},
{STONE, (120 / LAYER_DEPTH)},
{AIR, AIR_EXTERNAL}};
#define NUM_LAYERS (sizeof(layer)/sizeof(layer[0]))
const double del_t = (double) SECS_PER_DAY / TIME_INT_PER_DAY;
const double del_x = 1.0/(double) LAYERS_PER_M;
/*
* Compute an approximate steady-state initial T[*] vector based on a
* prorating of the INTERNAL_TEMP..EXTERNAL_TEMP range using r values.
* Returns the actual max count of the T[*].
*/
int initialSteadyStateApprox(double *T) {
double total_r, tmp_t;
int i, j, l;
/*
* Calculate the total R value for the wall profile.
*/
for (j=1, total_r= 0.0; jtype != AIR) {
/*
* In the case of conductive materials, the first layer is a
* material boundary, which is calculated as described in the
* intro, the other layers are internal so the standard W
* calculation is used.
*/
if (lp_1->type != AIR) { /* Material A | Material B boundary */
double k = 2.0 / (1.0/lp_1->K + 1.0/lp->K);
W[i] = (T[i-1] - T[i]) * k / del_x;
} else if (j==1) { /* Internal Air | Material boundary */
Tright_surface = (3*T[i] - T[i+1])/2;
W[i] = Q_FACTOR * (INTERNAL_TEMP - Tright_surface);
// Patch to inject two hours of a 2 deg temperature step.
if (l>639360 && l<=646560) W[i] +=2 * Q_FACTOR;
// if (W[i] > MAX_INTERNAL_HEAT) W[i] = MAX_INTERNAL_HEAT;
} else { /* Material A | Air gap | Material C boundary */
Tright_surface = (3*T[i] - T[i+1])/2;
Tleft_surface = (3*T[i-1] - T[i-2])/2;
W[i] = Q_FACTOR * (Tleft_surface - Tright_surface);
}
i++;
for (m=1; mK / del_x;
i++;
}
}
}
if (layer[NUM_LAYERS-1].type == AIR) { /* Material A | Air boundary */
Tleft_surface = (3*T[i-1] - T[i-2])/2;
W[i] = Q_FACTOR * (Tleft_surface - T[i]);
// Patch to inject two hours of sunshine
// if (l>639360 && l<=646560) W[i]-=400;
i++;
}
for (j=1, i=0; jtype != AIR) {
double q = del_t / (lp->rho * lp->cp * del_x);
for (m=0; m