# import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (6, 4)
import numpy as np
import pandas as pd
import fdsreader
from IPython.display import display
from PIL import Image
Data Analysis II#
This example demonstrates an analysis of slice data, here to determine the map of available safe egress time (ASET) and the temporal evolution of the smoke layer height. The used scenario is a multi-room appartment.
Fig. 3 Visualization of the appartment geometry.#
Full FDS input file
&HEAD CHID='Appartment' TITLE='Example in the Fire Simulation Lecture' /
&TIME T_END=300.0 /
&MISC BNDF_DEFAULT = F /
&BNDF QUANTITY = 'WALL TEMPERATURE' /
# ---------------------------------------
# Domain boundary
# ---------------------------------------
&VENT DB='YMAX', SURF_ID='OPEN' /
&VENT ID='OpenCeiling'
SURF_ID='OPEN'
XB=-1.5000,9.3000,4.3000,5.0000,2.2000,2.2000 /
&VENT ID='OpenLeft'
SURF_ID='OPEN'
XB=9.3000,9.3000,4.3000,5.0000,0.0000,2.2000 /
&VENT ID='OpenRight'
SURF_ID='OPEN'
XB=-1.5000,-1.5000,4.3000,5.0000,0.0000,2.2000 /
# ---------------------------------------
# Gas burner: methane
# ---------------------------------------
&SPEC ID = 'METHANE' /
&REAC FUEL = 'METHANE',
SOOT_YIELD = 0.060,
RADIATIVE_FRACTION = 0.02 /
# Burners:
&SURF ID = 'BURNER',
HRRPUA = 1000.0,
TMP_FRONT = 250.0,
COLOR = 'RASPBERRY' /
&OBST ID = 'BurnerBase',
SURF_IDS = 'BURNER', 'Wall Lining', 'Wall Lining',
XB = -0.20,0.20, -0.20,0.20, 0.00,0.20 /
# ---------------------------------------
# Surface definition (boundary conditions)
# ---------------------------------------
&SURF ID = 'Wall Lining',
COLOR = 'BEIGE',
DEFAULT = T,
THICKNESS(1:1) = 0.0254,
MATL_ID(1,1:1) = 'MARINITE',
MATL_MASS_FRACTION(1,1) = 1.0 /
&MATL ID = 'MARINITE'
EMISSIVITY = 0.9
DENSITY = 737.
SPECIFIC_HEAT = 1.2
CONDUCTIVITY = 0.12 /
--- Computational domain - 8 MESH
Cell Qty: 19008 | Size: 0.100·0.100·0.100m | Aspect: 1.0 | Poisson: No
&MESH ID='Domain' IJK=27,32,22 XB=-1.5000,1.2000,1.8000,5.0000,0.0000,2.2000 /
Cell Qty: 19008 | Size: 0.100·0.100·0.100m | Aspect: 1.0 | Poisson: No
&MESH ID='Domain.001' IJK=27,32,22 XB=1.2000,3.9000,1.8000,5.0000,0.0000,2.2000 /
Cell Qty: 19008 | Size: 0.100·0.100·0.100m | Aspect: 1.0 | Poisson: No
&MESH ID='Domain.002' IJK=27,32,22 XB=3.9000,6.6000,1.8000,5.0000,0.0000,2.2000 /
Cell Qty: 19008 | Size: 0.100·0.100·0.100m | Aspect: 1.0 | Poisson: No
&MESH ID='Domain.003' IJK=27,32,22 XB=6.6000,9.3000,1.8000,5.0000,0.0000,2.2000 /
Cell Qty: 19008 | Size: 0.100·0.100·0.100m | Aspect: 1.0 | Poisson: No
&MESH ID='Domain.004' IJK=27,32,22 XB=6.6000,9.3000,-1.4000,1.8000,0.0000,2.2000 /
Cell Qty: 19008 | Size: 0.100·0.100·0.100m | Aspect: 1.0 | Poisson: No
&MESH ID='Domain.005' IJK=27,32,22 XB=3.9000,6.6000,-1.4000,1.8000,0.0000,2.2000 /
Cell Qty: 19008 | Size: 0.100·0.100·0.100m | Aspect: 1.0 | Poisson: No
&MESH ID='Domain.006' IJK=27,32,22 XB=1.2000,3.9000,-1.4000,1.8000,0.0000,2.2000 /
Cell Qty: 19008 | Size: 0.100·0.100·0.100m | Aspect: 1.0 | Poisson: No
&MESH ID='Domain.007' IJK=27,32,22
XB=-1.5000,1.2000,-1.4000,1.8000,0.0000,2.2000 /
--- Geometric namelists from Blender Collections
-- Blender Collection: <Obstacles>
&OBST ID='AppartmentDoor',
XB=4.6000,5.1000,4.2000,4.3000,0.0000,2.2000 SURF_ID='Wall Lining' /
&OBST ID='AppartmentDoor.001',
XB=5.9000,6.4000,4.2000,4.3000,0.0000,2.2000 SURF_ID='Wall Lining' /
&OBST ID='AppartmentDoor.002',
XB=5.1000,5.9000,4.2000,4.3000,1.8000,2.2000 SURF_ID='Wall Lining' /
&OBST ID='Ceiling' SURF_ID='Wall Lining'
XB=-1.5000,4.6000,-1.5000,4.3000,2.2000,2.3100 SURF_ID='Wall Lining' /
&OBST ID='Ceiling.001' SURF_ID='Wall Lining'
XB=4.6000,9.4000,-1.5000,4.3000,2.2000,2.3100 SURF_ID='Wall Lining' /
&OBST ID='Floor', XB=-1.5000,9.4000,-1.5000,6.3100,-0.1100,0.0000
SURF_ID='Wall Lining' /
&OBST ID='WallBack',
XB=-1.5000,-1.4000,-1.4000,4.2000,0.0000,2.2000 SURF_ID='Wall Lining'
BNDF_OBST=T /
&OBST ID='WallBack.001',
XB=9.3000,9.4000,-1.5000,4.3000,0.0000,2.2000 SURF_ID='Wall Lining' /
&OBST ID='WallBathroom',
XB=6.5000,7.5000,1.4000,1.5100,0.0000,2.2000 SURF_ID='Wall Lining' /
&OBST ID='WallBathroom.001',
XB=8.3000,9.3200,1.4000,1.5100,0.0000,2.2000 SURF_ID='Wall Lining' /
&OBST ID='WallBathroom.002',
XB=7.5000,8.3000,1.4000,1.5100,1.8000,2.2000 SURF_ID='Wall Lining' /
&OBST ID='WallCorridorKitchen',
XB=6.4000,6.5000,0.4000,4.3000,0.0000,2.2000 SURF_ID='Wall Lining' /
&OBST ID='WallCorridorKitchen.001',
XB=6.4000,6.5000,-1.5000,-0.4000,0.0000,2.2000 SURF_ID='Wall Lining' /
&OBST ID='WallCorridorKitchen.002',
XB=6.4000,6.5000,-0.4000,0.4000,1.8000,2.2000 SURF_ID='Wall Lining' /
&OBST ID='WallCorridorLivingroom',
XB=4.5000,4.6000,-1.5000,-0.4000,0.0000,2.2000 SURF_ID='Wall Lining'
BNDF_OBST=T /
&OBST ID='WallCorridorLivingroom.001',
XB=4.5000,4.6000,0.4000,4.3000,0.0000,2.2000 SURF_ID='Wall Lining'
BNDF_OBST=T /
&OBST ID='WallCorridorLivingroom.002',
XB=4.5000,4.6000,-0.4000,0.4000,1.8000,2.2000 SURF_ID='Wall Lining'
BNDF_OBST=T /
&OBST ID='WallLeft',
XB=-1.5000,4.5000,-1.5000,-1.4000,0.0000,2.2000 SURF_ID='Wall Lining' /
&OBST ID='WallLeft.001',
XB=6.5000,9.3000,4.2000,4.3000,0.0000,2.2000 SURF_ID='Wall Lining' /
&OBST ID='WallLeft.002',
XB=4.6000,6.4000,-1.5000,-1.4000,0.0000,2.2000 SURF_ID='Wall Lining' /
&OBST ID='WallRight',
XB=-1.5000,4.5000,4.2000,4.3000,0.0000,2.2000 SURF_ID='Wall Lining' /
&OBST ID='WallRight.001',
XB=6.5000,9.3000,-1.5000,-1.4000,0.0000,2.2000 SURF_ID='Wall Lining' /
# ---------------------------------------
# Analytics
# ---------------------------------------
&SLCF ID = 'DoorLivingroomTempX',
PBX = 4.5500,
QUANTITY = 'TEMPERATURE',
CELL_CENTERED = T /
&SLCF ID = 'BurnerTempX',
PBX = 0.0000,
QUANTITY = 'TEMPERATURE'
CELL_CENTERED = T /
&SLCF ID = 'CorridorTempX'
PBX = 5.5000,
QUANTITY = 'TEMPERATURE',
CELL_CENTERED = T /
&SLCF ID = 'KitchenTempX'
PBX = 7.9000,
QUANTITY = 'TEMPERATURE',
CELL_CENTERED = T /
&SLCF ID = 'LivingRoomTempX'
PBX = 3.5000,
QUANTITY = 'TEMPERATURE',
CELL_CENTERED = T /
&SLCF ID = 'BurnerTempY'
PBY = 0.0000,
QUANTITY = 'TEMPERATURE',
CELL_CENTERED = T /
&SLCF ID = 'TempZ_1.0m'
PBZ = 1.0000,
QUANTITY = 'TEMPERATURE',
CELL_CENTERED = T /
&SLCF ID = 'TempZ_1.5m'
PBZ = 1.5000,
QUANTITY = 'TEMPERATURE',
CELL_CENTERED = T /
&SLCF ID = 'TempZ_2.0m'
PBZ = 2.0000,
QUANTITY = 'TEMPERATURE',
CELL_CENTERED = T /
&SLCF ID = 'DoorLivingroomSootDensityX',
PBX = 4.5500,
QUANTITY = 'DENSITY',
SPEC_ID = 'SOOT',
CELL_CENTERED = T /
&SLCF ID = 'BurnerSootDensityX',
PBX = 0.0000,
QUANTITY = 'DENSITY'
SPEC_ID = 'SOOT',
CELL_CENTERED = T /
&SLCF ID = 'CorridorSootDensityX'
PBX = 5.5000,
QUANTITY = 'DENSITY',
SPEC_ID = 'SOOT',
CELL_CENTERED = T /
&SLCF ID = 'KitchenSootDensityX'
PBX = 7.9000,
QUANTITY = 'DENSITY',
SPEC_ID = 'SOOT',
CELL_CENTERED = T /
&SLCF ID = 'LivingRoomSootDensityX'
PBX = 3.5000,
QUANTITY = 'DENSITY',
SPEC_ID = 'SOOT',
CELL_CENTERED = T /
&SLCF ID = 'BurnerSootDensityY'
PBY = 0.0000,
QUANTITY = 'DENSITY',
SPEC_ID = 'SOOT',
CELL_CENTERED = T /
&SLCF ID = 'SootDensityZ_1.0m'
PBZ = 1.0000,
QUANTITY = 'DENSITY',
SPEC_ID = 'SOOT',
CELL_CENTERED = T /
&SLCF ID = 'SootDensityZ_1.5m'
PBZ = 1.5000,
QUANTITY = 'DENSITY',
SPEC_ID = 'SOOT',
CELL_CENTERED = T /
&SLCF ID = 'SootDensityZ_2.0m'
PBZ = 2.0000,
QUANTITY = 'DENSITY',
SPEC_ID = 'SOOT',
CELL_CENTERED = T /
&SLCF ID = 'DoorLivingroomUVelX',
PBX = 4.5500,
QUANTITY = 'U-VELOCITY',
CELL_CENTERED = T /
&SLCF ID = 'DoorLivingroomDensityX',
PBX = 4.5500,
QUANTITY = 'DENSITY',
CELL_CENTERED = T /
&TAIL /
Precomputed simulation data
As shown in the introduction to fdsreader, the FDS simulation results can be read as Python data structures. The simulation consists out of eight meshes.
path_to_data = 'data/compartment/appartment_01/rundir'
sim = fdsreader.Simulation(path_to_data)
print(sim)
Simulation(chid=Appartment,
meshes=8,
obstructions=23,
slices=20,
smoke_3d=3)
ASET Map#
This section demonstrate how to compute ASET maps based on the procedure described in [SAS20]. The basic idea is to determine for each position in the compartment the first time that a tenability criterion is exceeded. The result is a time map.
Although the procedure is based on multiple criteria, for simlicity only one quantity, here the soot density, is used. Its values are evaluated in a z-normal slice at an arbitrary height of \(\SI{1.5}{\meter}\).
# get the soot density slice, normal to z at 1.5m height
slc = sim.slices.get_by_id('SootDensityZ_1.5m')
# as the simulation is based on multiple meshes, a global
# data structure is created, walls are represented as
# non-valid data points, i.e. nan
slc_data = slc.to_global(masked=True, fill=np.nan)
First, a visualisation of the data at a selected point is done with the imshow function.
# find the time index
it = slc.get_nearest_timestep(50)
# visualise the data
plt.imshow(slc_data[it,:,:].T, origin='lower', extent=slc.extent.as_list())
# add labels
plt.title(f'Soot Density at t={slc.times[it]:.2f}s')
plt.xlabel('position / m')
plt.ylabel('position / m')
plt.colorbar(orientation='horizontal', label=f'{slc.quantity.name} / {slc.quantity.unit}' )
# save output to file
plt.savefig('figs/appartment_soot_z.svg', bbox_inches='tight')
plt.close()
Fig. 4 Visualization of the soot density.#
Now, the local ASET values are computed:
Iterate over all spatial elements of the slice
Determine all points in time which exceed the tenability threshold
If this happens at any time, set the first time to be the local ASET value
# set arbitrary values as tenability threshold
soot_density_limit = 1e-4
# create a map with max ASET as default value
aset_map = np.full_like(slc_data[0], slc.times[-1])
# set walls to nan
aset_map[np.isnan(slc_data[0,:,:])] = np.nan
# 1D loop over all array indices, ix is a two dimensional index
for ix in np.ndindex(aset_map.shape):
# find spatialy local values which exceed the given limit
local_aset = np.where(slc_data[:, ix[0], ix[1]] > soot_density_limit)[0]
# if any value exists
if len(local_aset) > 0:
# use the first, i.e. first in time, as the local ASET value
aset_map[ix] = slc.times[local_aset[0]]
With the computed map, a graphical respresentation of the ASET map is done the same way as with the other quantities. Here, a discrete color map is used.
# create a discrete (12 values) color map
# cmap = matplotlib.cm.get_cmap('jet_r', 12)
cmap = plt.cm.get_cmap('jet_r', 12)
# visualise the data
plt.imshow(aset_map.T, origin='lower', extent=slc.extent.as_list(), cmap=cmap)
plt.title(f'ASET Map with Soot Density Limit of {soot_density_limit:.1e}')
plt.xlabel('x position / m')
plt.ylabel('y position / m')
plt.colorbar(orientation='horizontal', label='time / s' );
# save output to file
plt.savefig('figs/appartment_aset_map.svg', bbox_inches='tight')
plt.close()
Fig. 5 ASET map for the outlined scenario.#
Smoke layer#
In this example, the smoke layer height is analysed. The distinction made here is based on a simple threshold in temperature: The local smoke layer height is given by the lowest point above a given temperature. The evaluation is done based on a slice across the burner and normal to the x-direction.
# find the slice
slc = sim.slices.get_by_id('BurnerTempX')
# convert it to a global data structure and get the coordinates
slc_data, slc_coords = slc.to_global(masked=True, fill=np.nan, return_coordinates=True)
First, the data at a arbitrary point in time is visualsied. The white parts represent the obsticles.
# pick a time index
it = slc.get_nearest_timestep(150)
# visualise the data
plt.imshow(slc_data[it,:,:].T, origin='lower', vmax=200, extent=slc.extent.as_list())
plt.title(f'Temperature at t={slc.times[it]:.2f}s')
plt.xlabel('y position / m')
plt.ylabel('z position / m')
plt.colorbar(orientation='horizontal', label=f'{slc.quantity.name} / {slc.quantity.unit}' )
# save output to file
plt.savefig('figs/appartment_temp_slice.svg', bbox_inches='tight')
plt.close()
Fig. 6 Temperature slice across the burner.#
Now, for each y-position the z-indices are found, where the temperature exceedes the limit temperature. The lowest value is the smoke layer height at the y-position.
# set temperature limit
temperature_limit = 75
# create a data array to store the local height values, default
# is the maximal z-coordinate
layer_height = np.full(slc_data.shape[1], slc_coords['z'][-1])
# loop over all indices
for ix in range(len(layer_height)):
# find indices which exceed the limit
lt = np.where(slc_data[it, ix, :] > temperature_limit)[0]
# if there are any, pick the lowest one
if len(lt) > 0:
layer_height[ix] = slc_coords['z'][lt[0]]
The resulting values can now be plotted over the slice file, to check for plausibility.
# slice data
plt.imshow(slc_data[it,:,:].T, origin='lower', vmax=200, extent=slc.extent.as_list())
plt.title(f'Temperature at t={slc.times[it]:.2f}s')
plt.xlabel('y position / m')
plt.ylabel('z position / m')
plt.colorbar(orientation='horizontal', label=f'{slc.quantity.name} / {slc.quantity.unit}' );
# smoke layer height
plt.plot(slc_coords['y'], layer_height, '.-', color='red')
# save output to file
plt.savefig('figs/appartment_temp_slice_height.svg', bbox_inches='tight')
plt.close()
Fig. 7 Temperature slice with the local smoke layer height.#
Using the above approach for a single point in time, a loop over all times can be used to compute, e.g., the mean and standard deviation of the smoke layer height.
layer_mean = np.zeros_like(slc.times)
layer_stddev = np.zeros_like(slc.times)
res = np.zeros(slc_data.shape[1])
for it in range(len(slc.times)):
res[:] = slc_coords['z'][-1]
for ix in range(len(res)):
lt = np.where(slc_data[it, ix, :] > temperature_limit)[0]
if len(lt) > 0:
res[ix] = slc_coords['z'][lt[0]]
layer_mean[it] = np.mean(res)
layer_stddev[it] = np.std(res)
# plot the mean and stddev values as functions of time
plt.plot(slc.times, layer_mean, label='Mean Smoke Layer Height')
plt.plot(slc.times, layer_stddev, label='Stddev of Smoke Layer Height')
plt.grid()
plt.legend()
plt.xlabel('Time / s')
plt.ylabel('Height / m')
# save output to file
plt.savefig('figs/appartment_layer_mean_stddev.svg', bbox_inches='tight')
plt.close()
Text(0, 0.5, 'Height / m')
Both values can be combined and visualised jointly, where the standard deviation is used to indicate a fluctuation band around the mean value.
# plot the mean
plt.plot(slc.times, layer_mean, label='Mean Smoke Layer Height')
# plot a band around the mean, using the stddev as band borders
plt.fill_between(slc.times, layer_mean-layer_stddev, layer_mean+layer_stddev, color='C0', alpha=0.3)
# show the floor for reference
plt.ylim(bottom=0)
plt.grid()
plt.legend()
plt.xlabel('Time / s')
plt.ylabel('Height / m')
# save output to file
plt.savefig('figs/appartment_layer_mean_band.svg', bbox_inches='tight')
plt.close()
Fig. 8 Temporal development of the mean and stdandard deviation of the smoke layer height.#
If parts of the region shall be excluded in the analysis, a coordinate dependent mask can be used for this.
# find indices, where the y coordinate is between the given values
ymin = 1
ymax = 4
coord_mask = np.where((slc_coords['y'] > ymin) & (slc_coords['y'] < ymax))
# slice data
plt.imshow(slc_data[it,:,:].T, origin='lower', vmax=200, extent=slc.extent.as_list())
plt.title(f'Temperature at t={slc.times[it]:.2f}s')
plt.xlabel('y position / m')
plt.ylabel('z position / m')
plt.colorbar(orientation='horizontal', label=f'{slc.quantity.name} / {slc.quantity.unit}' );
# smoke layer height
plt.plot(slc_coords['y'][coord_mask], layer_height[coord_mask], '.-', color='red')
# save output to file
plt.savefig('figs/appartment_temp_slice_height_mask.svg', bbox_inches='tight')
plt.close()
Fig. 9 Temperature slice with the local smoke layer height, limited to the masked region.#
The above procedure can be reused, yet the computation of the mean and standard deviation is carried out on the masked values.
for it in range(len(slc.times)):
res[:] = slc_coords['z'][-1]
for ix in np.ndindex(res.shape):
lt = np.where(slc_data[it, ix, :] > temperature_limit)[1]
if len(lt) > 0:
res[ix] = slc_coords['z'][lt[0]]
# computation is carried out on the masked values now
layer_mean[it] = np.mean(res[coord_mask])
layer_stddev[it] = np.std(res[coord_mask])
# same plot as above
plt.plot(slc.times, layer_mean, label='Mean Smoke Layer Height')
plt.fill_between(slc.times, layer_mean-layer_stddev, layer_mean+layer_stddev, color='C0', alpha=0.3)
plt.ylim(bottom=0)
plt.grid()
plt.legend()
plt.xlabel('Time / s')
plt.ylabel('Height / m')
# save output to file
plt.savefig('figs/appartment_layer_mean_band_mask.svg', bbox_inches='tight')
plt.close()
Fig. 10 Temporal development of the mean and stdandard deviation of the smoke layer height in the masked region.#