Multiphysics Modelling of Stone Wool Fire Resistance

Abstract

In fire resistance tests, stone wool’s organic matter undergoes exothermic oxidative reactions sustained by external heat, causing mass transfer in the structure. The previous modelling attempts, lacking the mass transfer physics, fall short in predicting the temperature of high density and high organic content samples. To fill this gap in the fire engineering modelling capability, we include mass transfer in our calculation, and validate the model using experimental fire resistance data. As an alternative, we use a heat conduction -based model lacking the gas transfer but with reaction kinetics coupled to the stone wool’s organic mass %. The results show that the thermal effects of the oxidative degradation can be predicted by introducing the simplified diffusion processes. The oxygen transfer and exothermic reactions depend upon the amount of organic content, and the uncertainty of temperature predictions is $\pm \,20\%$. In average, temperatures and critical times are more accurately predicted by the heat conduction model, while, the peak temperature prediction uncertainty is low ($\pm \,10\%$) with the multiphysics model. The uncertainty compensation method reduces the difference between the two model predictions. Nevertheless, further validation study is needed to generalize the uncertainty compensation metrics. Finally, we demonstrate how a gas flow barrier on the cold side (sandwich) can effectively reduce the peak temperature of the high organic content-stone wools.

1 Introduction

Stone wools are widely used in residential buildings and industrial facilities for thermal insulation and fire resistance [1]. Under conventional safety regulation, the design of such passive fire protection is based on the standard fire test carried out according to EN 1363-1. The test procedure for thermal insulation performance consists of measuring the cold-side temperature of the sample exposed to a combustion chamber depicting the ISO-834 standard time-temperature curve [2].

The stone wool thermal response against the fire is a multiphysics problem. Being a porous structure, the heat transfer occurs in all three modes: conduction, convection and radiation [3]. Also, the wool usually contains organic material that undergoes chemical degradation. The process is slow and flameless, similar to smouldering combustion [4], but sustained by the furnace heating rather than the reaction itself. The chemical process triggers mass transfer within the structure [5], and as a consequence, the availability of oxygen within the porous region affects the rate of the chemical reactions and the overall heating behaviour of the structure [6].

To support and extend the stone wool thermal resistance testing, there have been several modelling attempts. Lee [7] investigated the effect of fibre orientation on the intra-fibre radiative heat transfer, observing that the radiative fluxes are highest when the fibres are perpendicular to the boundary. Andersen and Dyrboel [8] predicted the radiative fluxes using Planck mean properties and compared them to the one obtained using spectral and flux-weighted properties. They showed that the Planck mean property-based predictions are accurate for wavelengths higher than $10\,\upmu\hbox {m}$ and highly porous stone wools. Veiseh and Hakkaki-Fard [9] estimated the thermal conductivity by numerically solving the steady-state radiative and conductive heat transfer, and compared it to the experimentally measured values. They found that the predicted conductivities excellently resembled the measured one for low-density stone wools.

Some studies have included the effect of organic material degradation on the overall heat transfer. For example, Karamanos et al. [10] formulated a theoretical approach of updating thermal conductivity as a function of temperature and organic material concentration for different chemical and structural configurations of the stone wools. They showed that the wool properties are permanently changed upon exceeding a certain temperature threshold and that the insulation degrades when the temperature is increased even further. Livkiss et al. [11] and Andres et al. [6] carried out various tests and studied the chemical degradation and thermal conductance of stone wools. They found that the thermogravimetric analysis (TGA) and micro-scale combustion calorimetry (MCC) tests show three distinct mass loss rate peaks, the first two indicating the degradation of organic matter and the last one that of fibres. They then formulated the finite-difference solutions to 1D heat transfer and chemical reactions. The solutions, however, were accurate only for low density and low organic content stone wools, indicating that the TGA and MCC tests in a simple way could not be used to model the reaction heat of high density and high organic content stone wools. Also, suggesting that the modelling of oxygen controlled reaction scheme could explain the stone wool heating behaviour at a larger scale.

The previous modelling attempts have reproduced the termination of exothermic combustion due to the low temperature and material consumption, but have not explicitly taken into account the oxygen availability. Two possible mechanisms behind the oxygen controlled degradation can be identified: the limitation of transport due to the low porosity (high density) and oxygen replacement by the high yield of gaseous oxidation products (high organic content). In this paper, we adopt the hypothesis of oxygen controlled degradation rate and test it by formulating a 1D multiphysics model of the heat and mass transfer, coupled with the chemical decomposition, and by solving it numerically for a range of stone wool properties. As a reference model, representing the classical approach, we will use a heat conduction -based model which lacks porous media gas transfer but with optimized decomposition kinetics. We will compare the modelling results against experimental data to validate the models and to conclude if the oxygen control is, indeed, the limiting factor for the decomposition of the organic content.

Due to the uncertainties in the fire conditions and stone wools material properties, the fire performance evaluation needs to be probabilistic [12, 13]. We estimate the probabilistic cold-side temperatures considering the fire condition and stone wool material as the input stochastic. We also study the model sensitivities to different material and boundary conditions.

2 Method

2.1 Fire Resistance Test

Figure 1 depicts small scale fire resistance test. A sample of size $60 \,\hbox{cm}\times 60\,\hbox {cm}$ having thickness l is separated from the combustion chamber using 1 mm thin steel plate. The cold side is exposed to the ambient ($T_\infty$) and the hot side follows the ISO 834 standard fire curve,

$$\begin{aligned} T_{{\text {h}}}(t) = T_\infty + 345{{\text {log}}}_{10}(8t + 1). \end{aligned}$$

(1)

Figure 1

Test geometry for small scale fire resistance test, EN 1363-1. a Sample size and thermocouple locations, b cross-section, c cold-side temperature measurements

Table 1 Material and geometrical properties of the investigated stone wools

Altogether, 30 tests (one sample per material) were performed with varying wool density, thickness and organic content. Table 1 lists the material and geometrical properties of the investigated stone wool types. The dimensions are according to EN 12085 standard. The fourth column shows the Loss On Ignition (LOI) values according to EN 13820. It represents the organic content of stone wool (as % weight). The fifth column shows the intact fibre % measured using a method adapted from ASTM C1335-12. It represents the amount of fibrous material of the total inorganic solids (% weight), assumed to impact the thermal insulation properties via thermal radiation. The sixth and seventh columns show fibre mean thickness, $d_f$, and orientation, $\theta$, measured using the internal analysis methods of material manufacturer. The last column shows the Rosseland mean extinction coefficient ($\beta$), calculated using intact fibre %, mean thickness ($d_f$) and orientation ($\theta$) [10],

$$\begin{aligned} \beta&=\rho e, \\ e&=\left( 133.0-11.23d_f+0.5d_f^2-0.0024d_f^3\right) \cdot \frac{{{\text {cos}}}(\theta )}{2}, \nonumber \\ \rho&=(\rho _{{{\text {w}}},0}-\rho _{{{\text {o}}},0})\cdot ({{\text {intact\, fibre }}} \%), \nonumber \end{aligned}$$

(2)

where $\rho$ is the material density, e is the specific extinction coefficient, and $\rho _{{{\text {w}}},0}$ and $\rho _{{{\text {o,s}}}}$ are initial density for stone wool and organic matter respectively.

The test installations are according to the European standard EN 1363-1. The furnace was gas-fuelled, and the temperatures were recorded using K-type thermocouples. Figure 1c shows the cold side installation of the thermocouples. The measuring ends of the thermocouples are covered with $30 \times 30 \times 2 \,\hbox {mm}$ inorganic insulating pads. The pad and thermocouple were attached to the slab by heat-resistant glue or pins. Figure 2 shows the recording of furnace temperature and cold side surface temperatures of all the tests. The furnace temperature meets the allowed tolerance specified in EN 1363-1. The cold side temperatures are the averages of 8 thermocouples.

Figure 2

Test measurements. Left: furnace temperature. Right: cold-side temperatures (the average of 8 thermocouples)

2.2 Multiphysics Model

The stone wool with density $\rho _{{{\text {w}}}}$ consists of fibres ($\rho _{{{\text {f}}}}$), organic material ($\rho _{{{\text {o}}}}$), products of the decomposed organic matter ($\rho _{{{\text {p}}}}$), and air ($\rho _{{{\text {a}}}}$),

$$\begin{aligned} \rho _{{{\text {w}}}}=\rho _{{{\text {f}}}}+\rho _{{{\text {o}}}} +\rho _{{{\text {a}}}}+\rho _{{{\text {p}}}}. \end{aligned}$$

(3)

Figure 3

Multiphysics model computational set up. $l_s$ and l are the steel plate and stone wool thicknesses, respectively. The gap in the stone wool indicates that the lengths are not in scale ($l>>l_s$)

The mass transfer in the stone wool domain triggers due to the chemical decomposition of the organic matter. The direct contribution of the decomposition on the overall energy balance is negligible, but significant heat may be released if the decomposition product oxidice at elevated temperature [5]. As the information of the exact chemical composition and reaction steps are not available, we assume an approximate single step scheme consisting of two lumped gas species, air and products,

$$\begin{aligned} {{\text {C}}}_8{{\text {H}}}_6{{\text {O}}}_2 + 5.5(\overbrace{{{\text {O}}}_2+3.76{{\text {N}}}_2} ^{\begin{array}{c} {{\text {Air}}} \end{array}})\rightarrow \overbrace{ 5{{\text {CO}}}_2+20.7{{\text {N}}}_2+3{{\text {H}}}_2{{\text {O}}}}^{\begin{array}{c} {{\text {Products}}} \end{array}}+3{{\text {C}}}, \end{aligned}$$

(4)

where 1g of organic matter (${{\text {C}}}_8{{\text {H}}}_6{{\text {O}}}_2$, phenol-formaldehyde resin) reacting with 5.63 g air produces 6.45 g gas products and 0.27 g char.

The mass transfer is assumed to be purely diffusive. The solved governing equations and their initial conditions, therefore, are the continuity equation for organic mass, diffusion equation for the products mass, and the diffusion equation for the enthalpy in the slab:

$$\begin{aligned} \frac{\partial \rho _{{\text {o}}}}{\partial t}&=-{\dot{\omega }}\rho _{{{\text {w}}}},\quad 0\le x\le l,\quad \rho _{{{\text {o}}},0} ={{\text {LOI}}}\cdot \rho _{{{\text {w}}}}, \end{aligned}$$

(5)

$$\begin{aligned} \frac{\partial \rho _{{\text {p}}}}{\partial t}&=\frac{\partial }{\partial x}\left( D\frac{\partial \rho _{{\text {p}}}}{\partial x}\right) +{\dot{\omega }}\rho _{{{\text {w}}}}f,\quad 0\le x\le l, \quad \rho _{{{\text {p}}},0} =0, \end{aligned}$$

(6)

$$\begin{aligned} \rho c\frac{\partial T}{\partial t}&=\frac{\partial }{\partial x}\left( k\frac{\partial T}{\partial x}\right) +{\dot{\omega }}\rho _{{{\text {w}}}}\Delta h,\quad -l_s\le x\le l, \quad T_0 =T_\infty . \end{aligned}$$

(7)

The factor ($f=6.45$) is according to the decomposition scheme Eq. 4. The heat of combustion ($\Delta h$) is 25 MJ/kg for the organic material according to an EN ISO 1716 standard measurement. The rate of the oxidation reaction ${\dot{\omega }}$ is assumed to follow an Arrhenius form [14]

$$\begin{aligned} {\dot{\omega }}=Ae^{\left( -\frac{E_{{{\text {a}}}}}{{{\text {R}}}T}\right) } \left( \frac{\rho _{{{\text {o}}}}}{\rho _{{{\text {w}}}}} \right) ^nX_{{{\text {O}}}_2}^{n_{{{\text {O}}}_2}}. \end{aligned}$$

(8)

Here, $X_{{{\text {O}}}_2}$ is the oxygen volume fraction calculated using the volume proportions of the air and products,

$$\begin{aligned} X_{{{\text {O}}}_2}=0.21X_{{\text {a}}},\quad X_{{\text {a}}}=1-X_{{\text {p}}},\quad X_{{\text {p}}}=\frac{\rho _{{\text {p}}} \rho _{{{\text {a,s}}}}}{\rho _{{{\text {p,s}}}}}, \end{aligned}$$

(9)

where 0.21 is the oxygen volume fraction in air. $\rho _{{{\text {a,s}}}}$ and $\rho _{{{\text {p,s}}}}$ are the specific density of air and products respectively, see “Appendix A”. Figure 3 shows the domain specific implementation of the governing equations. The mass conservation is considered only for the stone wool domain, whereas the energy conservation is solved for the entire slab, consisting of the steel plate (thickness $l_s= 1 \,\hbox {mm}$) and wool (thickness l).

Regarding the boundary conditions, we assume that the gas diffusion is blocked on the hot-side, depends upon the surface diffusion coefficient ($D_2$) on the cold side, and the radiative and convective heat fluxes affect on both sides,

$$\begin{aligned} -D\frac{\partial \rho _{{\text {p}}}}{\partial x}&= 0 \quad {{\text { at }}}\, x =0, \end{aligned}$$

(10)

$$\begin{aligned} D\frac{\partial \rho _{{\text {p}}}}{\partial x}&= -D_2\rho _{{\text {p}}} \quad {{\text { at }}}\, x =l, \end{aligned}$$

(11)

$$\begin{aligned} -k\frac{\partial T}{\partial x}&= h_1(T_{{{\text {h}}}} - T(x)) + \epsilon _\Gamma \sigma (T_{{{\text {h}}}}^4 - T(x)^4) \quad {{\text { at }}}\, x =-l_s, \end{aligned}$$

(12)

$$\begin{aligned} -k\frac{\partial T}{\partial x}&= h_2(T(x)-T_{\infty }) + \epsilon _\Gamma \sigma (T(x)^4 - T_{\infty }^4) \quad {{\text { at }}}\, x =l. \end{aligned}$$

(13)

The generic material properties of steel, air, and the products are presented in Appendix. The stone wool bulk properties depend on T or $\rho _{{{\text {o}}}}$, or both, and are calculated by assuming the medium between the fibres to be air. This is because the change in the proportions of air and the products negligibly affects the bulk properties. Similarly, the effect of char is assumed to be negligible because on average, the char at the end of tests is less than 1.5%. The bulk thermal conductivity is according to Karamanos approach [10], in which, the thermal conductivity of the pure fibre is a temperature dependent expression having coefficients $a_0,\ldots a_3$ [15]. The equivalent expression includes conductive and radiative part [16]. The conductive part includes the parallel, || , and the perpendicular, $\perp$, components, that depend upon porosity ($\phi$). The radiative part is calculated using the refractive index value $n=1.2$, Stefan–Boltzmann constant ($\sigma$), and extinction coefficient ($\beta$).

$$\begin{aligned} \begin{array}{lr} \begin{aligned} k_{{\text {w}}}&{}= k_{{\text {w}}}(T,\rho _{{{\text {o}}}},\theta ,\beta ),\\ k_{{\text {f}}}&{}=a_0+a_1T^{-2}-a_2T^{\frac{1}{2}}+a_3T\\ \phi &{}= 1-\frac{\rho _{{{\text {o}}}}}{\rho _{{{\text {o,s}}}}} -\frac{\rho _{{{\text {w}}}}-\rho _{{{\text {o}}}}}{\rho _{{{\text {f,s}}}}}. \\ k_{||}&{}=\phi k_{{\text {a}}}+(1-\phi )k_{{\text {f}}}, \\ k_{\perp }&{}=\frac{k_{{\text {a}}}k_{{\text {f}}}}{{\phi } k_{{\text {f}}}+(1-\phi )k_{{\text {a}}}},\\ k_{{\text {cond}}}&{}=k_{||}\sin (\theta )^2+k_{\perp }\cos (\theta )^2,\\ k_{{\text {radi}}}&{}=\frac{16n^2\sigma T^3}{3\beta },\\ k_{{\text {w}}}&{}=k_{{\text {cond}}}+k_{{\text {radi}}}. \end{aligned}&{}\quad \begin{array}{lr} \hbox {Constants} &{} \hbox {values} \\ a_0&{}1.98{{\text {E}}}{+}0 \\ a_1&{}7.26{{\text {E}}}{+}3 \\ a_2&{}6.45{{\text {E}}}{-}2 \\ a_3&{}1.27{{\text {E}}}{-}3 \\ \rho _{{{\text {o,s}}}}&{}1300 \hbox {kg~m}^{-3} \\ \rho _{{{\text {f,s}}}}&{}2850 \hbox {kg~m}^{-3}\\ \end{array} \end{array} \end{aligned}$$

(14)

The bulk heat capacity is the weighted fraction of fiber, organic matter, and air, in which, the heat capacity of the pure fiber is a temperature-dependent expression having coefficients $a_0,\ldots ,a_2$ [15]. The equivalent expression uses the density of the solids ($\rho _{{\text {s}}}$) obtained by subtracting the density of the gas,

$$\begin{aligned} \begin{array}{ll} \begin{aligned} c_{{\text {w}}}&{}= c_{{\text {w}}}(T,\rho _{{{\text {o}}}}), \\ c_{{\text {f}}}&{}=a_0+a_1T-a_2T^{-2}, \\ \rho _{{\text {s}}}&{}=\rho _{{{\text {w}}}}-\phi \rho _{{{\text {a,s}}}},\\ c_{{\text {w}}}&{}=\frac{\rho _{{\text {s}}}}{\rho _{{{\text {w}}}}^2} \left[ (\rho _{{{\text {w}}}}-\rho _{{{\text {o}}}})\cdot c_{{\text {f}}} +\rho _{{{\text {o}}}}\cdot c_{{\text {o}}}\right] +\frac{1}{\rho _{{{\text {w}}}}}(\rho _{{{\text {w}}}} -\rho _{{\text {s}}})\cdot c_{{\text {a}}}. \end{aligned} &{}\quad \begin{array}{lr} \hbox {Constants} &{} \hbox {values} \\ a_0&{}9.50{{\text {E}}}{+}2 \\ a_1&{}1.60{{\text {E}}}{-}1 \\ a_2&{}1.35{{\text {E}}}{-}1 \\ c_{{\text {o}}}&{}0.2 \,\hbox {Jkg}^{-1}\hbox {K}^{-1} \\ \end{array} \end{array} \end{aligned}$$

(15)

where, $c_{{\text {a}}}$ is the heat capacity of air, see “Appendix A”. Figure 4 shows $k_{{\text {w}}}$ and $c_{{\text {w}}}$ variation to temperature. The other dependencies are constant to the initial state of the stone wool, See Table 1.

Figure 4

Temperature dependence of thermal conductivity (left) and heat capacity (right)

The heat transfer coefficient on the cold side ($h_2$) is qualitatively close to the natural convection of a vertical wall, as shown in the left plot Fig. 5. On the hot-side($h_1$), Wang [17] suggests value between 2 Wm⁻² K⁻¹ and 60 Wm⁻² K⁻¹. Using this range we carry out random test and opt $h_1=18.5\,\hbox {Wm}^{-2}~\hbox {K}^{-1}$ for minimum error in the cold side temperature. The gas diffusivity (D) is calculated using Chapman and Enskog equation, and shown in the right plot in Fig. 5 [18]. The gas diffusion coefficient on the cold-side ($D_2$) is assumed to be $10^{-3}$ order of $h_2$.

Figure 5

Left: cold-side heat transfer coefficient. Right: average diffusivity of the gaseous mixture, D [18]

The chemical decomposition parameters (A, $E_{{\text {a}}}$, n and $n_{{{\text {o}}}_2}$) are optimized by carrying out a Monte-Carlo (MC) simulation of $\hbox {N}=1000$ samples using Latin Hypercube Sampling (LHS) technique, and finding a minimum error in the cold-side temperatures. The sample range is listed in Table 2 and the search criterion is the mean square error. The third column shows the optimized values.

Table 2 The sample spaces and the optimum values for the chemical decomposition parameters

The numerical implementation is using the Finite Element Method (FEM) and a custom Python script. The finite elements are linear and uniform of size 1 mm. The time discretization follows the explicit Euler method. The time step is adaptive to the maximum temperature rise of 1 °C. The solver adopts a time-step between a millisecond to about a hundred seconds. The temperature and density-dependent material matrices are updated every step to ensure the correct field-dependent properties.

The numerical convergence was tested by comparing the results with reduced element size (0.5 mm). The maximum change in the cold-side temperature was less than 0.3%. The predictions were also benchmarked using COMSOL Multiphysics. The maximum difference in temperature prediction was $\simeq 3\%$. The computational cost with a custom Python script is high. The execution time on intel Core i7 CPU (quadcore) with 2.73 GHz, was nearly one hour for the simulation of one hour of stone wool heating. The same scenario in the COMSOL Multiphysics took about a minute. The custom Python script, however, was less burdensome for distributed computing, hence, was used in the MC simulation that was carried out for the optimization of reaction parameters.

2.3 Heat Conduction Model

A separate model was constructed using Fire Dynamics Simulator (FDS, version 6.7.0), a commonly used tool for fire-driven flows and pyrolysis of structures [19]. The model solves parabolic heat diffusion equation coupled to the Arrhenius equation,

$$\begin{aligned} \rho c\frac{\partial T(x,t)}{\partial t}= \frac{\partial }{\partial x}\left( k\frac{\partial T(x,t)}{\partial x}\right) +\rho _{{{\text {w}}}}{\dot{\omega }}\Delta h,\quad {\dot{\omega }}={{\text {A}}}e^{\left( -\frac{{{\text {E}}}_{{{\text {a}}}}}{{{\text {R}}}T} \right) }\left( \frac{\rho _{{{\text {o}}}}}{\rho _{{{\text {w}}}}}\right) ^{n}. \end{aligned}$$

(16)

The boundary conditions are same as Eq. 13. The stone wool bulk properties are defined to change only to temperature, i.e., the same as Fig. 4.

Figure 6

Steps for the optimization of the chemical decomposition parameters: A, $E_{{\text {a}}}$ and n

Figure 7

Optimization of chemical decomposition parameters. The dots represents optimum values for each wool types. The dashed lines represent the optimized function or a value for all wool types

Table 3 Optimized values of the chemical decomposition parameters for all the stone wool types

The chemical decomposition parameters (A, $E_{{\text {a}}}$ and n) are correlated to stone wools’ organic material mass %. The correlation functions are formulated based on the pattern seen in test data. Figure 6 shows the steps for the formulation of the correlation functions. First, MC simulations were carried out for N=1000 samples of A, $E_{{\text {a}}}$ and n, generated using the LHS technique and the range listed in Table 2. Then, comparing the least square error in the cold-side temperatures, the optimal A, $E_{{\text {a}}}$ and n values were selected for each wool types. The correlations were, then, searched by plotting the optimized values (A, $E_{{\text {a}}}$ and n) against LOI, density and thickness. No correlation was found between A, $E_{{\text {a}}}$ and n and the stone wool density or thickness. Similarly, no correlation was found between $E_\text {a}$ and stone wool properties. The only correlation found was between A and LOI, and n and LOI. The left and right plot in Fig. 7 show the correlations, where A and n are the second-order polynomial and logarithmic expansion of LOI, respectively. Then the second round MC simulation was carried out to optimize the coefficients of the correlation functions, $a_0, a_1, a_2, b_0, b_1$, and $E_{{\text {a}}}$, where $a_0, a_1, a_2$ and $b_0, b_1$ are the coefficients for A and n respectively. The shaded region in the plots of Fig. 7 represents the searched spaces. The dotted lines are the optimized functions, with parameters listed in Table 3.

2.4 Sensitivity and Stochastic Analysis

2.4.1 Sensitivity to Wool Properties and Boundary Condition

The cold-side temperature sensitivity to stone wool properties was studied for both types of models using three test pairs, each having a difference in either thickness, density or LOI value, see Table 4. The pairs were selected in such a way that there is a significant difference in the investigated property and negligible differences in the remaining two properties.

Table 4 The pairs of test cases to investigate the cold-side temperature sensitivity to stone wool properties

Figure 8

The time-temperature curve representing the three different hot-side boundary conditions

Similarly, the sensitivity to hot side temperature history was studied with three different hot side boundary conditions (Fig. 8): The ISO 834 standard curve, Eq. 1, constant value, $T_{{\text {g}}} =945 \,^\circ \hbox {C}$, i.e., same as the ISO 834 standard curve at 60 min, and an exponential curve, $T_{{\text {g}}}=T_\infty +\hbox {exp}(\hbox {t}/8.65)$, i.e., the mirror of the ISO 834 standard curve.

Additionally, to study the effect of oxygen transfer, the cold-side temperatures were simulated for a different ratio of gas diffusion and heat transfer coefficient, $\frac{D_2}{h_2}=2{{\text {E}}}-4,6{{\text {E}}}-4,\ldots ,1.8{{\text {E}}}-3$. The significance of the oxygen transfer was further confirmed by carrying out two additional experiments with stone wool containing 4.8% organic matter. In the first test, the steel plate was placed only on the hot side, as before. In the second test, steel plates were placed on both sides.

2.4.2 Stochastic Analysis and Uncertainty Management

To form a general view of the stone wool behaviour and model performance within a wide range of engineering applications, we estimate the wool cold side temperature for a randomly chosen stone wool material and a range of time-temperature curves depicting various compartment fire scenarios. Figure 9 shows the time-temperature curves, the location of the stone wool layer, and the geometry of the compartment. The floor and ceiling are concrete, and the walls are Gypsum board, except one corner where the wall is made of stone wool. Table 5 lists the material properties. The curves represent office building compartment fire according to Annex A of EN 1991-1-2:2000.

Figure 9

Upper left: a parametric time-temperature curve representing various compartment fire scenarios. Upper right: location of the stone wool protective layer. Lower: a fire compartment with enclosure size: $1.6a\times a\times 3.0 \,{{\text {m}}}^3$ and opening size: $1.2a\times h \,{{\text {m}}}^2$

Table 5 Material properties and thickness of the compartment boundaries in the stochastic analysis

Table 6 The input parameter distributions for the fire scenarios of the stochastic analysis

The random inputs are the room width, a, height of the opening, h, fire load density, $q_{{\text {f,t}}}$, the thermal diffusivity of walls, $\alpha _{{\text {w}}}$, and the thermal diffusivity of ceiling or floor, $\alpha _{{\text {f}}}$. The samples size is $N =1000$, and the sampling method is LHS. Table 6 lists their mean and standard deviation or the lower and upper values. The last variable in the table, i.e., Wool type, is to select a type of stone wool material from 30 different varieties listed in Table 1.

Figure 10

Schematic diagram showing the steps for stochastic analysis and uncertainty management

From the viewpoint of the probabilistic applications of the simulation models, we are interested in knowing if the two model prediction match after the uncertainty correction. Figure 10 shows the steps for the uncertainty correction. First, the model uncertainties are quantified in terms of the systematic bias, $\delta$, and the second central moment of the random errors, ${\widetilde{\sigma }}_\epsilon$. The former one defines the average variation in the predicted, ${\hat{\psi }}_i$, and measured, $\psi _i$, quantities, and the latter one defines the average variation in the difference of the predicted and measured quantities, for M number of validation cases [20, 21],

$$\begin{aligned} \delta = \frac{\mu _{{\hat{\psi }}}}{\mu _\psi },\quad {{\text {and}}}\quad {\widetilde{\sigma }}_\epsilon =\frac{1}{\mu _{{\hat{\psi }}}} \left[ \frac{1}{M-1}\sum _{i=1}^{M}\left( {\hat{\psi }}_i-\delta \cdot \psi _i\right) ^2\right] ^\frac{1}{2}. \end{aligned}$$

(17)

Next, we apply the previously developed uncertainty compensation model [21], and $\delta$ and ${\widetilde{\sigma }}_{\epsilon }$, to calculate corrected outputs,

$$\begin{aligned} \psi _i = \frac{1}{\delta } \left[ \mu _{{\hat{\psi }}} + \left( {\hat{\psi }}_i-\mu _{{\hat{\psi }}}\right) \sqrt{1-\left( \frac{\mu _{{\hat{\psi }}}{\widetilde{\sigma }}_{\epsilon }}{\sigma _{{\hat{\psi }}}}\right) ^2}\right] ,\quad i=1,2, .., N, \end{aligned}$$

(18)

where ${\hat{\psi }}_i$ is the original simulation output, and $\mu _{{\hat{\psi }}}$ and $\sigma _{{\hat{\psi }}}$ are the mean and standard deviation (all i).

3 Results

3.1 Validation

3.1.1 Cold-side temperature

The fire resistance tests were simulated and the cold-side temperatures were predicted using both multiphysics and heat conduction model. Figure 11 compares the measured and predicted cold-side temperatures for wool 7, 22 and 30, whose heating behaviour differs from each-other: wool 7 has low organic content, wool 22 and 30 have high organic content, and wool 22 differ from 30 in terms of density (see Table 1). The bell-shaped area in the plots represents the region affected by the reaction heat. The prediction in this region is comparatively more accurate for the low-density stone wool, i.e., wool 30. For the low organic content, i.e., wool 7, the prediction is more accurate with the heat conduction model. For high density and high organic content, i.e., wool 22, the prediction is more accurate with the multiphysics model. On average, the thermal effect of reaction heat is captured well by the multiphysics model.

Figure 11

The evolution of measured and predicted cold-side temperatures

Figure 12

The scatter plots of the measured and predicted cold-side temperatures

Figure 12 shows the scatter plots of the peak values and the temperatures at 30 and 60 minutes obtained from all the tests. The scatter values are aligned along the diagonal for both model types, and the models have similar performance in the prediction of the temperature at a specific time. For the peak temperatures, however, the multiphysics model predictions are slightly closer to the diagonal.

3.1.2 Time to critical temperature

To evaluate the models’ capability to predict fire resistance times, the times of the cold side rising above an arbitrary critical temperature $T_{{\text {Cr}}}$ were determined from the time-temperature curves. The plots in Fig. 13 show these times as functions of the normalized critical temperature, having value one at $T_\infty + 140\,^\circ \hbox {C}$, which is the highest acceptable mean surface temperature in EN 1363-1. The plots show that the times are underpredicted by the multiphysics model and overpredicted by the heat conduction model. For low organic content, i.e., wool 7, the heat conduction model predictions are more accurate. For high organic content, i.e., wool 22 and 30, the multiphysics model predictions are more accurate. There is an unexpected delay in the rising of wool 30 cold-side temperature predicted by the heat transfer model. This leads to a higher discrepancy in the predicted times.

Figure 13

The measured and predicted times at which the cold side temperature rises to $T_{{\text {Cr}}}\,^\circ \hbox {C}$

Figure 14

The scatter plots of the times at which the cold-side rises by 100 °C, 140 °C and 180 °C

Figure 14 shows the scatter plots of the times at which the cold side temperature rises by 100 °C, 140 °C and 180 °C. The right-hand-side plots do not show the result of all the tests because their cold side temperatures did not reach the 180 °C threshold. The plots show that the values for the heat conduction and multiphysics models are slightly above and below the diagonal, respectively. The heat conduction model predictions are closer to the diagonal.

3.1.3 Modelling Uncertainty

Figure 15 shows the $\delta$ and ${\widetilde{\sigma }}_\epsilon$ for the temperature and time predictions. $\hbox {N}=30$ for all the listed outputs except “time to $T_\infty +140^\circ \hbox {C}$” and “$T_\infty +180^\circ \hbox {C}$”. For these quantities, $\hbox {N}=28$ and $\hbox {N}=8$ respectively, see Fig. 14. The last two are the average values of the outputs. The plot indicates that the multiphysics model over-predicts temperatures and under-predicts times, while the heat conduction model over-predicts times.

Figure 15

Modelling uncertainty, $\delta \pm {\widetilde{\sigma }}_\epsilon$, in temperature and time prediction

3.2 Qualitative analysis

Using the multiphysics model, the field variables (temperature, organic mass %, and oxygen availability) were estimated along the thickness of the stone wools. The plots in Fig. 16 show the profiles of these quantities for wools 7, 22 and 30, at times 10, 20 and 40 min. The temperature plots (first row) show how the cold side temperatures start to grow after 10 min. At 20 min, the cold side temperature gradients of the wools 7 and 22 become similar to the hot side gradients and exceed the hot size gradients for wool 30. Low-density-wools have a higher thermal conductivity at high temperatures than the high-density wool, and thus the heat is transferred easier.

Figure 16

Temperatures, organic content mass %, and oxygen availability along with the thickness of wool 7, 22 and 30, and at 10, 20 and 40 min

The second row shows the organic matter mass concentration. The plots indicate that the oxidative degradation propagates from the hotter-side to the colder-side. For wool 22, some portion of the organic matter remains unreacted on the hotter side because there is no oxygen available (lower middle plot) for oxidative decomposition, due to the high generation of product gas (high organic content sample). For wools 7 and 30, the organic matter decomposition is uninterrupted. Some degree of oxygen limitation is observed in these wools as well (lower, left and right plot) but that does not last for long because the product gas generation is low (low organic content sample).

The third row shows the oxygen availability. The y-axis value 1 means that the porous media is fully air, and 0 means that the air is not available. The plots indicate that for the wool 7 (left plot), the product gas never covers the entire porous region, i.e., the oxygen is always available. The time-wise change in the curves is most significant for wool 30 (middle and right plots). This means the faster air recovery for a wool with comparatively smaller organic mass %.

3.3 Sensitivity and Stochastic Analysis

3.3.1 Sensitivity to Material Properties

Figure 17 shows the cold side temperature sensitivity to the change in stone wool material properties. The plots in the first two rows indicate that the shift of the time-temperature curve in time towards higher times is due to the increase in either thickness or the density of the stone wool. The higher value of thickness or density means a longer delay in the rising of the cold-side temperature. Such delay is predicted by both types of model. Nevertheless, the ones predicted by the multiphysics model are closer to the measured ones.

Figure 17

Cold-side temperature sensitivity to stone wool thickness (upper row), density (middle row) and the LOI value (lower row)

The plots in the third row of Fig. 17 indicate that the increase in LOI value increases the peak of the time-temperature curve. Such an effect is predicted by both types of models. With the heat conduction model, however, there is an unexpected shift in time towards the right side. This might be because of improper $k_{{\text {w}}}$ modelling. When the stone wool heats up, the organic mass % decreases and contributes to the increase of $k_{{\text {w}}}$, see Eq. 14. In the heat conduction model, $k_{{\text {w}}}$ only depends on the temperature. As the thermal conductivity is not increasing with the organic mass % decrease, the temperature curve is shifted later.

3.3.2 Sensitivity to Boundary Condition

Figure 18 shows the cold side temperature sensitivity to different hot side temperature history. The results indicate that the constant boundary condition amplifies the models’ temperature response. It also decreases the discrepancy between two types of models. With the exponential hot-side boundary, the stone wools temperature remains low, but the model predictions are significantly different in an absolute sense.

Figure 18

Cold-side temperature sensitivity to different boundary conditions: Iso standard curve (upper row), a constant value (middle row), exponential curve (lower row)

Figure 19 shows the cold-side temperatures for different ratio of gas diffusion and heat transfer coefficient, $\frac{D_2}{h_2}=2{{\text {E}}}-4,6{{\text {E}}}-4,\ldots ,1.8{{\text {E}}}-3$. The plots indicate that the temperatures are insensitive to the gas diffusion for wool 7. This is because of the low organic content. The curve peak, however, decreases with an increase in the gas diffusion for wool 22 and 30. This means that, for the high organic content stone wool, the high fire resistivity (low-temperature peak) can be achieved by reducing the gas transfer.

Figure 19

Cold-side temperature sensitivity to the ratio of gas diffusion and heat transfer coefficient, $\frac{D_2}{h_2}$

Figure 20

Sensitivity to open or closed boundary. Left: Cold-side temperatures. Middle: organic mass % at the end of test. Right: visual inspection of the decomposition zone at the end of test

Figure 20 shows the measured and predicted (multiphysics model) cold side temperatures for the open and closed cold side boundary. The temperature peak is not observed for the wool with closed boundaries. In both bases, the predicted temperatures are mostly in good agreement with the measured temperatures. However, after the temperature peak in the open-boundary case, the nearly steady-state temperature is slightly under-predicted by the model. The middle plot in Fig. 20 shows the predicted organic content at the end of the tests. It indicates that for the open boundary wool, the decomposition occurs throughout the thickness, but reaches only up to $\sim \,20$ mm for the closed one. This prediction resembles the decomposition zone seen in the photograph of the samples taken after dissecting them into half. The photograph shows also that the thickness of the open boundary -wool (left sample) has decreased more than that of the closed one (right sample).

3.3.3 Stochastic Analysis

Figure 21 shows the predicted cold-side temperature solution spaces for both stone wool models. The dotted lines represents the $\Phi =0.1$, 0.5 and 0.9 fractile values of the distribution. The plots indicate that the bell-shaped region, in the case of the multiphysics model, covers a relatively wide time range ($\sim \,10$ to $\sim \,40 \,\hbox {min}$). With the heat conduction model, the oxidative reactions start later, but after the initiation, the reactions take place faster, covering a narrower range of time, and ending by $\sim 30\,\hbox {min}$.

Figure 21

The cold-side temperature solution space obtained from the stochastic analysis. The dotted curves represent $\phi =0.1$, 0.5 and 0.9 fractile value of the distribution

Figure 22

The cumulative density functions for temperature and time prediction. The dotted and continuous line respectively represent the multiphysics and heat conduction models

Figure 23

The corrected cumulative density functions for temperature and time prediction. The dotted and continuous line respectively represent the multiphysics and heat conduction models

Figure 22 compares the Cumulative Density Functions (CDF) for six different output variables: peak temperature, temperatures at 30 min and 60 min, and times for three different threshold temperatures corresponding to the temperature increase by 100 °C, 140 °C, and 180 °C. The time CDF curves for the 140 and 180 degree temperature increase end below unity because these threshold temperatures were reached in only a fraction of the simulated fires. These distributions were calculated by first collecting the statistics of the cases where the criterion was met, and then normalizing the curve with the share of those cases.

The plots show that there can be large discrepancies in the time distributions despite the small discrepancy in the peak temperature distributions. For most output quantities, the differences between the model distributions comply with the test uncertainty metrics in Figure 15. The only exception is “Time to $T_\infty +180^\circ \hbox {C}$”, for which the discrepancy in the distributions is larger than the discrepancy in the biases shown in Fig. 15. The reason is probably related to the small size of the data behind the validation metrics.

Figure 23 shows the distribution after uncertainty correction. Based on the visual evaluation, we can say that the best improvement is seen in “Time to $T_\infty +140^\circ \hbox {C}$”. The results indicate that the uncertainty correction is most effective for the outputs which are quantitatively but not qualitatively different, and for which the difference of the predicted statistics complies with the the validation metrics (Fig. 15).

4 Discussion

We reported the modelling uncertainties of heat conduction and multiphysics models, showing that the overall prediction is more accurate with the heat conduction model. However, the heat conduction model predictions would be far less accurate without correlating the chemical degradation parameters to the LOI values of the stone wools. Also, note here that the model sensitivities to the changes in stone wool material properties and hot-side boundary conditions are more realistic with the multiphysics model. In real applications, the choice of the modelling method also depends on the computational cost. In this respect, the alternative of the LOI-correlated decomposition scheme in the heat conduction method is more efficient than to solve the entire mass diffusion problem.

The stochastic analysis with a range of time-temperature curves and the full range of different model types showed noticeable differences between the two models, Fig. 22. The discrepancies reduced (Fig. 23) after utilizing the previously developed error compensation method, Eq. 18. The compensation, however, is not perfect because the correction parameters (Fig. 15) do not fully match with the discrepancies seen in the stochastic outputs (Fig. 22). In other words, the validation cases are not generic enough: The stone wool properties (Table 1) do not vary systematically, and the fire condition is identical in all the tests. For a highly non-linear problem as the heat transfer of porous, reactive material with many dependent variables (organic mass %, density, thickness, fibre orientation, intact fibre %), the 30 validation tests were too few for the generalization of the model uncertainty. New tests with systematic variations in material and fire conditions would be needed to generalize the uncertainties.

The temperature curves obtained for the open or closed boundary systems, Fig. 20, clearly demonstrated the requirement of the external oxygen flow for the high exothermic temperature peaks. This finding can be used for analyzing the benefits of sandwich -type structures over open-surface wool products. The results also revealed that the predicted temperatures in the end of the experiment are lower than the measured ones, and the difference is greater for the sample with open cold side boundary. The visual inspection of the photograph showed the reduction in sample thickness, which is not taken into account in the models. The cause for the thickness reduction is the fibre material decomposition or the structural deformation of the fibre. In fact, by numerical experiments we found out that the model could accurately predict the observed final temperature if the sample thickness reduction was explicitly specified in the model. The future modelling efforts may consider including the mechanical response in order to improve the temperature predictions.

The factor that determines the oxygen unavailability in the stone wool is the amount of organic matter. Livkiss et al. [11] observed that ignoring the dependency of degradation reactions on the oxygen availability leads to inaccurate temperature estimation for high-density and high-organic content stone wools. They concluded that the low porosity is the cause of the oxygen unavailability in the high-density stone wools. Considering the wools studied in this work, the stone wool porosity is high ($\ge 0.95$) regardless of the stone wool density (38 to $147\, \hbox {kg~m}^{-3}$). The high density, thus, does not necessarily affect the oxygen transfer. One example is Wool 7 ($100\, \hbox {kg~m}^{-3}$), for which the temperature profiles are unchanged for different values of gas diffusion coefficients, see Figure 19. The confusion arises because the amount of organic matter may be the same despite the differences in either density or LOI. For example, the amount of organic matter in $100\, \hbox {kg~m}^{-3}$ stone wool with 5% LOI is same as in $50\, \hbox {kg~m}^{-3}$ wool with 10% LOI. It seems that the lack of oxygen transfer is seen either in the samples with high LOI or the high-density samples with sufficient LOI.

We assume a purely diffusive mass transfer and explain the heating behaviour of high organic content and high-density stone wools. The validity of the purely diffusive mass transfer-assumption can be examined by calculating the ratio of the convective and diffusive time scales [22] $\tau = (\mu L^2 \Delta p^{-1} B_0^{-1})/(L^2 D^{-1})= {\mu D}(\Delta p B_0)^{-1}$, where the permeability coefficient $B_0 \sim 10^{-10}\,\hbox {m}^2$ [23], dynamic viscosity $\mu \sim 10^{-5}\,\hbox {Pa}\cdot \hbox {s}$, and gas diffusivity $D\sim 10^{-5}\,\hbox {m}^2\hbox {s}^{-1}$. With pressure differences $\Delta p> 10\,\hbox {Pa}$, $\tau < 1$, i.e. the convection dominates over diffusion. However, as the convection of expanded gas and combustion products away from the reaction zone and the diffusion of oxygen from the cold-side surface to the reaction zone act in opposite directions, we assumed that the diffusion process alone can yield first-order predictions of the oxygen availability. The positive results support this assumption, but the convective flow should be considered in future research.

5 Conclusion

In this study, we presented a multiphysics approach of modelling stone wool heating behaviour, and an alternative heat conduction -based approach lacking the gas transfer but coupled to the LOI-dependent reaction kinetics. The capability of the models to predict the fire resistance thermal insulation of different stone wools was studied by validating against experimental data and by sensitivity and stochastic analyses.

The multiphysics model shows that the stone wool temperature depends upon the availability of oxygen. In stone wools with high organic content, the oxygen entrainment, the chemical decomposition and the release of heat may be interrupted due to the high rate of product gas production, whereas, in low organic content stone wools, this is not the case. The high density of stone wools, however, cannot be claimed to be the cause for porosity-related oxygen unavailability because all the investigated stone wools, despite variation in the density, are highly porous ($\ge 0.95$) and open for oxygen transfer.

Both types of models were found to be capable of reproducing experimentally observed bell-shaped, temperature peaks in the case of high organic content-stone wools and the continuous temperature profiles otherwise, with $\pm \, 20\%$ uncertainty. In average, the temperature and critical times prediction are more accurate with the heat conduction model, while, the peak temperature prediction uncertainty is low ($\pm \, 10\%$) with the multiphysics model. Note that the heat conduction model performs poorly without the LOI -correlated reaction parameters.

Similarly, both models predicted the expected responses to variations in material composition and hot-side boundary conditions. The constant boundary condition amplifies the stone-wool response, the change in thickness or the density of the stone wool shifts the time-temperature curve in time, and the LOI value determines the time-temperature curve peak of high organic content -stone wools. The accurate prediction of the shifting of the time-temperature curve requires tracking of the change of stone wool organic material mass % and updating of the continuum model conductivity and capacity matrices at appropriate time intervals.

The models’ stochastic performances are notably different near the bell-shaped temperature region, where the multiphysics models predicts wider temperature distribution than the heat conduction model. Correcting the models using the uncertainty metric from the validation set reduces the discrepancies. Nevertheless, for perfect correction, the uncertainty metric would need further generalization with validation tests covering a wider range of material properties and fire conditions.

Considering the use of open and sandwich -type stone wool products in fire protection, we observed that the peak cold side temperature of the high organic content -stone wool could be efectively reduced by preventing the flow on the cold side. A thin layer of non-permeable material, thus, can be used on the cold side to block oxygen transfer and to increase the fire resistivity of high organic content stone wools.

Appendix: Generic material properties

$$\begin{aligned} \rho _{{\text {steel}}}&=7850 [{{\text { kg/m}}}^{3}]\nonumber \\ k_{{\text {steel}}}(T)&=6.4E1-3.5E-2T [{{\text {W}}}/({{\text {m}}}\cdot {{\text {K}}})]\;{{\text { for }}}\; T{{\text { in \,Kelvin}}}\nonumber \\ c_{{\text {steel}}}(T)&=4.3E2+7.7E-1T-1.69E-3T^2+2.22E-6T^3 [{{\text {J/Kg}}}/^\circ {{\text {C}}}]\;{{\text { for }}}\; T<600^\circ {{\text {C}}} \nonumber \\&=6.66E2+1.3E4/(738-T)[{{\text {J/Kg}}}/^\circ {{\text {C}}}]\;{{\text { for }}}\; T>=600.0 {{\text { and }}} T<735^\circ {{\text {C}}} \nonumber \\&=5.5E2+1.8E4/(T-731.0)[{{\text {J/Kg}}}/^\circ {{\text {C}}}]\;{{\text { for }}}\; T>=735 \;{{\text { and }}}\; T<900^\circ {{\text {C}}}\nonumber \\&=650.0 [{{\text {J/Kg}}}/^\circ {{\text {C}}}] ,\quad \;{{\text {for }}}\; T>900^\circ {{\text {C}}}\\ \rho _{{{\text {a,s}}}}&=350.56T^ {-0.999} [{{\text {J/Kg}}}/^\circ {{\text {C}}}]\;{{\text { for }}}\; T{{\text { in \,Kelvin}}}\\ k_{{\text {a}}}(T)&=4.3E-3\left[ c_1T_r^{-1}+c_2T_r^{-2.0/3.0}+ c_3T_r^{-1.0/3.0}+c_4+c_5T_r^{1.0/3.0}+c_6T_r^{2.0/3.0}+c_7T_r+c_8T_r^{4.0/3.0}\right. \\&\quad \left. +\,c_9T_r^{5.0/3.0}\right] [{{\text {W}}}/({{\text {m}}}\cdot {{\text {K}}})],\quad c1=33.97,\quad c2=-164.70,\quad c3=262.11,\quad c4=-21.53 \\&\quad c5=-443.45,\quad c6=607.34,\quad c7=-368.80,\quad c8=111.30,\quad c9=-13.41,\; {{\text{and}}}\\&\quad T_r=T/132.52\;{{\text { for }}}\; T\; {{\text{in \, Kelvin}}} \\ c_{{\text {a}}}(T)&=-2E-7T^3-4E-4T^2-1.45E-1T+1.01E3[{{\text {J/Kg}}}/{{\text {K}}}]\;{{\text { for }}}\; T\; {{\text {in\, Kelvin}}}\\ \rho _{{{\text {p,s}}}}&=\frac{5\rho _{{{\text {CO}}}_2{{\text {,s}}}}+20.7\rho _{{{\text {NO}}}_2{{\text {,s}}}}+3\rho _{{{\text {H}}}_2{{\text {O,s}}}}}{5+20.7+3}\\ \rho _{{{\text {CO}}}_2{{\text {,s}}}}&=1.98[{{\text { kg/m}}}^{3}]\\ \rho _{{{\text {NO}}}_2{{\text {,s}}}}&=1.25[{{\text { kg/m}}}^{3}]\\ \rho _{{{\text {H}}}_2{{\text {O,s}}}}&=0.60[{{\text { kg/m}}}^{3}]\\ X_{O_2}&=0.21\left( \phi -\rho _{{{\text {p}}}}\frac{\rho _{{{\text {a,s}}}}}{\rho _{{{\text {p,s}}}}}\right) =0.21\left( \phi -\rho _{{{\text {p}}}}\frac{1.225}{1.31}\right) \end{aligned}$$