Finite Element Analysis is the technique of dividing a large object into much smaller pieces or elements, modelling those small elements and then recombining the model results to give a result for the whole. The beauty of the approach is that while modelling the whole object may be difficult or numerically intractable, the elemental level model can be a simplification of the real world phenomenon and yet still yield good results when the elements are recombined.
In the example below we look at heat loss in a steel column. In domestic properties, steel columns and beams are often used to provide structural support. In comparison to traditional materials such as wood and stone, steel has a high strength to volume ratio. The modern I beam crossection is strong in compression and in shear. This allows the creation of large open living areas that just would not be possible otherwise. Steel is however a good conductor of heat. The question is therefore, when there is pressure to reduce domestic energy consumption for heating, what price is being paid for the luxury of open plan living, and what can be done about it.
The two dimensional model considered here, is for a steel column of 3.5 metres height. It is assumed that the edges upper 3 metres are pinned at the internal room temperature of 20C and the lower 0.5metres at an external temperature of 0C. Note, this means we are assuming the heating system will deliver *whatever it takes* to maintain the interior at 20C. Also note there are no conduction or radiation edge effects and the floor is assumed to be of negligible thickness.
Our starting point is the well known heat equation. This says that the rate of change of temperature at any point is proportional to the sum of the rates of change of the rate of change in temperature in the x direction and y direction i.e.
du/dt = alpha * (d2u/dx + d2u/dy)
In the above, alpha, the proportionality constant is the thermal conductivity of the material (measured in W/mK) divided by the product of the material density (Kg/m3) and heat capacity (J/KgK). The equation quantifies what is easy to understand intuitively - a difference in temperature between an object and its surroundings leads to an inflow of energy which raises an object's temperature in proportion to the object's heat capacity. Heat capacity is best understood by considering how it varies among common materials: metals have relatively low heat capacities - they heat up quickly when a heat source is applied (e.g. frying pan). It requires 480J of energy to raise the temperature of a block of steel by one degree. Brick and stone by contrast require double the energy to achieve the same temperature increase. Water has one of the highest common heat capacities 4200J/kgK (which is why it is such an excellent material for energy transport in heating and cooling systems).
We will represent the steel as a 350x50 Numpy matrix (17,500 elements). Numpy is chosen for its speed of calculation - the vectorisation of many matrix operations means our programme will run much faster than if we were to use native Python data objects. Additionally the mathematical orientation of the language helps express the problem in code in a succint functional manner. The benefits of succint, expressive, easy to understand code, can not be overstated.
Lets set up our environment, we'll imagine a 4m x 1m modelling arena within which our column will sit
import numpy as np
real_length=4.0
real_width=1.0
length = 400
width = 100
dx = real_length / length
dy = real_width / width
internal_temp = 20
external_temp = 0
steel = {'density': 7800, 'conductivity':50, 'heat_capacity': 480}
Now lets set up some arrays to model the properties of interest. We'll use booleans to indicate where a material is and is not and also where edges sit (edges will become important later when we want to calculate the net heat flows into and out of the steel).
u = np.full((length, width), internal_temp, dtype=float)
density = np.full_like(u, 0)
conductivity = np.full_like(u, 0)
heat_capacity = np.full_like(u, 0)
steel_mask = np.full_like(u, False, dtype=bool)
exterior_mask = np.full_like(u, False, dtype=bool)
interior_mask = np.full_like(u, False, dtype=bool)
interior_xplus_edge = np.full_like(u, False, dtype=bool)
interior_xminus_edge = np.full_like(u, False, dtype=bool)
exterior_xplus_edge = np.full_like(u, False, dtype=bool)
exterior_xminus_edge = np.full_like(u, False, dtype=bool)
interior_yplus_edge = np.full_like(u, False, dtype=bool)
interior_yminus_edge = np.full_like(u, False, dtype=bool)
exterior_yplus_edge = np.full_like(u, False, dtype=bool)
exterior_yminus_edge = np.full_like(u, False, dtype=bool)
setting out our shapes, namely the steel, its edges, the room interior (at 20C) and exterior (at 0C)
interior_mask[0:320,:] = True
exterior_mask[320:,:] = True
steel_mask[20:370,20:70] = True
interior_xplus_edge[20:320, 20]=True
interior_xminus_edge[20:320, 69]=True
interior_yplus_edge[20,20:70] = True
exterior_xplus_edge[320:370, 20] = True
exterior_xminus_edge[320:370, 69] = True
exterior_yplus_edge[369,20:70] = True
interior_mask = np.logical_and(interior_mask, ~(np.logical_or(steel_mask)))
exterior_mask = np.logical_and(exterior_mask, ~(np.logical_or(steel_mask)))
density[steel_mask] = steel['density']
conductivity[steel_mask] = steel['conductivity']
heat_capacity[steel_mask] = steel['heat_capacity']
density[insulation_mask] = insulation['density']
conductivity[insulation_mask] = insulation['conductivity']
heat_capacity[insulation_mask]= insulation['heat_capacity']
heat_capacity[exterior_mask] = 9999
density[exterior_mask] = 9999
conductivity[exterior_mask] = 9999
heat_capacity[interior_mask] = 9999
density[interior_mask] = 9999
conductivity[interior_mask] = 9999
In discretising the differential terms of the heat equation we can approximate the gradients in temperature in the x and y directions so that for example du/dx = (u[x + dx] - u[x]) / dx. To get d2u/dx2 we repeat the process using du/dx rather than u[x]. With some light rearranging this yields
d2y/dx2 = (u[x+1,y] - u[x,y])/dx**2 + (u[x-1,y] - u[x-1,y])/dx**2
and similarly
d2u/dy2 = (u[x,y+1] - u[x,y])/dy**2 + (u[x,y-1] - u[x,y])/dy**2
On the other side of the heat equation, discretization gives us
du/dt = (u[t + dt,x,y] - u[t,x,y])/dt
Putting this all together and rearranging gives
u[t+dt,x y] = u[t,x,y]) + dt * (
alpha * (u[x+1,y] - u[x,y])/dx**2
+ alpha * (u[x-1,y] - u[x-1,y])/dx**2
+ alpha * (u[x,y+1] - u[x,y])/dy**2
+ alpha * (u[x,y-1] - u[x,y])/dy**2 )
Note that this is a recurrence relation - if we have the temperature for a particular point at a particular time u[t,x,y] we can use the formula to find the temperature a moment later, u[t+dt, x,y]. This will form the basis of our simulation.
Three modelling issues to consider:
The main iteration loop is then
vol_heat_capacity = density * heat_capacity
alpha = conductivity / vol_heat_capacity #units Wm2/kJ
max_dt = np.nanmin(1/(1*alpha*(1/dx**2 + 1/dy**2)))
dt = min(1, max_dt)
for i in np.arange(100000):
#set boundary conditions - assume edges are pinned at temperature
# and accept whatever is the implied heat flow.
#3.5m steel column, upper 3m is within heated envelope,
#lower 0.5m is outside heated envelope.
#Aassume underfloor provides no shelter from external conditions.
u[interior_mask] = internal_temp
u[exterior_mask] = external_temp
#we have a heterogenous collection of materials therefore the thermal conductivity
#relevant to each boundary is a combination of the conductivities on either side
#of each boundary (not sure - use minimum or mean?)
conductivity_left = np.minimum(conductivity[2:, 1:-1] , conductivity[1:-1, 1:-1])
conductivity_right = np.minimum(conductivity[:-2, 1:-1] , conductivity[1:-1, 1:-1])
conductivity_bottom = np.minimum(conductivity[1:-1, 2: ] , conductivity[1:-1, 1:-1])
conductivity_top = np.minimum(conductivity[1:-1, :-2] , conductivity[1:-1, 1:-1])
#temperature increase is the sum of increases from the four cardinal directions
#tempered by the heat capacities of the materials
u[1:-1, 1:-1] = u[1:-1, 1:-1] + dt * (
(u[2:, 1:-1] - u[1:-1, 1:-1]) * conductivity_left / (vol_heat_capacity[1:-1, 1:-1] * dx**2)
+ (u[:-2, 1:-1] - u[1:-1, 1:-1]) * conductivity_right / (vol_heat_capacity[1:-1, 1:-1] * dx**2)
+ (u[1:-1, 2: ] - u[1:-1, 1:-1]) * conductivity_bottom/ (vol_heat_capacity[1:-1, 1:-1] * dy**2)
+ (u[1:-1, :-2] - u[1:-1, 1:-1]) * conductivity_top / (vol_heat_capacity[1:-1, 1:-1] * dy**2)
)
We'll use matplotlib to display the temparature distribution
import matplotlib.pyplot as plt
plt.pcolormesh(u[::-1], cmap=plt.cm.jet, vmin=u.min(), vmax=u.max())
plt.gca().add_line (line = plt.Line2D((20, 20), ( 30,380), lw=1, color='grey'))
plt.gca().add_line (line = plt.Line2D((20, 70), ( 30, 30), lw=1, color='grey'))
plt.gca().add_line (line = plt.Line2D((20, 70), (380,380), lw=1, color='grey'))
plt.gca().add_line (line = plt.Line2D((70, 70), ( 30,380), lw=1, color='grey'))
plt.colorbar()
The high thermal conductivity of the steel means that the interior portion stabilises at internal temperature, the exterior at the external temperature with a relatively small transition zone in between.
gradient_x_plus = -np.minimum(conductivity[:,1:], conductivity[:,:-1]) * (u[:,:-1] - u[:, 1:])
gradient_x_minus = -np.minimum(conductivity[:,1:], conductivity[:,:-1]) * (u[:,1: ] - u[:, :-1])
gradient_y_plus = -np.minimum(conductivity[1:,:], conductivity[:-1,:]) * (u[:-1,:] - u[1:, :])
gradient_y_minus = -np.minimum(conductivity[1:,:], conductivity[:-1,:]) * (u[1:,: ] - u[:-1, :])
internal_heat_loss = (np.sum(np.abs(gradient_x_plus[interior_xplus_edge[:,1:]]))
+ np.sum(np.abs(gradient_x_minus[interior_xminus_edge[:,:-1]]))
+ np.sum(np.abs(gradient_y_plus[interior_yplus_edge[1:,:]]))
+ np.sum(np.abs(gradient_y_minus[interior_yminus_edge[:-1,:]]))
)
external_heat_loss = (np.sum(np.abs(gradient_x_plus[exterior_xplus_edge[:,1:]]))
+ np.sum(np.abs(gradient_x_minus[exterior_xminus_edge[:,:-1]]))
+ np.sum(np.abs(gradient_y_plus[exterior_yplus_edge[:-1,:]]))
+ np.sum(np.abs(gradient_y_minus[exterior_yminus_edge[1:,:]]))
)
Running the code gives the answer, 2700W heatloss. This is per metre thickness in the unmodelled third dimension. If we say roughly that a typical steel will have flanges and web of 10mm, then the total thickness is 30mm or 0.03metres, so the expected heatloss is 2700x0.03 = 81W. For a room with four such columns that's 325W. To put this into context, a typical double glazed window might have a heat loss of 1W/m2K and a single glazed window 2W/m2K. Our steels are the equivalent of 16m2 of double glazing. So by poor design we have erased most of the benefit of say converting a house to double glazing. What is also a concern, is that if we look at the temperature of the column at the floor junction (inspect the matrix u) it is 6C. If the rest of the room is at 50% with relative humidity at 20C, then a quick check of a psychrometric chart shows the air will saturate at just under 10C, so at 6C we will see condensation forming, with all the concomittent issues of mold growth, fabric damage and health risks for the occupants. So what can be done. Let's investigate two possible solutions.
insulation = {'density': 115, 'conductivity':0.05, 'heat_capacity': 1000}
We'll modify our arrays to take account of 100mm insulation added
insulation_mask = np.full_like(u, False, dtype=bool)
insulation_mask[370:380,20:70] = True
insulation_mask[320:380,10:20] = True
insulation_mask[320:380,70:80] = True
interior_xplus_edge[20:320, 20] = True
interior_xminus_edge[20:320, 69] = True
interior_yplus_edge[20,20:70] = True
interior_yplus_edge[320,10:20] = True
interior_yplus_edge[320,70:80] = True
exterior_yplus_edge[379,10:80] = True
exterior_xplus_edge[320:380,10] = True
exterior_xminus_edge[320:380,79] = True
interior_mask = np.logical_and(interior_mask, ~(np.logical_or(steel_mask, insulation_mask)))
exterior_mask = np.logical_and(exterior_mask, ~(np.logical_or(steel_mask, insulation_mask)))
conductivity[insulation_mask] = insulation['conductivity']
heat_capacity[insulation_mask] = insulation['heat_capacity']
From the picture we can immediately understand the benefit of insulation. Heat loss is reduced from 2700W/m to 17W/m. The problem floor temperature rises to a fraction under 20C.
However such structural insulation would be expensive, perhaps prohibitively so, not to mention the problem of trying to find skills in the domestic sector for work normally only employed on commercial scale projects. Therefore as an alternative, let's look at the effect of insulating the sides of the column only. This would be a simpler solution - there is no complex engineering and the materials and skills are readily available.
Heat loss is 156W/m (or 5W for our four columns of 10mm flange/web thickness). Not quite as good as the rolls-royce solution in absolute terms, but pretty impressive ontheless and a much better return on investment.
Note that we have not wrung all we could out of the simulation results. We have only looked at the end state. Since, we have accurately modelled the heat flows and material properties, we could stop the iteration at any point to look at rate of spread of heat through the material. As the drive for economy pushes levels of insulation up, and heating power down, it is possible to reach a situation where changes in state driven by timer controls and external weather conditions mean a comfortable steady state is never achieved [much to the chagrin of the occupant who might have spent considerable sums on say the latest heat pump technology]. Analysing heat flow profiles over time to understand when this is a risk is a sensible precursor to any capital spend. For an example analysis showing how small a heating system can be before the heat-up wait time becomes problematic - see here
The problem of discretising a differential equation and solving via iteration is one that comes up in a wide variety of disciplines - physics, chemistry, engineering, finance. See here for an example from finance.
Back to other example analyses.