An Application of Relative Entropy in Structural Safety Analysis of Elastoplastic Beam under Fire Conditions

Abstract

The main aim of this work is to propose a new algorithm of reliability assessment for steel civil engineering structures subjected to fire temperatures. This new algorithm is based upon the relative probabilistic entropy concept elaborated by Bhattacharyya, and this probabilistic distance is sought in-between extreme and admissible deformations of some structural beam subjected to higher temperatures. Similar to the First Order Reliability Method, this approach uses the first two probabilistic characteristics of the structural response, when structural output may be modelled with the use of Gaussian distribution. The probabilistic structural response is found here using hybrid computational technique–the Finite Element Method system ABAQUS with its fully coupled thermo-elastic analysis with 3D solid elements and probabilistic modules implemented in the computer algebra system MAPLE. The probabilistic response is determined via a triple stochastic analysis, which is based on the classical Monte-Carlo simulation, iterative generalized stochastic perturbation technique, and also using semi-analytical probabilistic calculus. Final determination of the relative entropy in the Serviceability Limit State of the given structure and its comparison with the results obtained using the FORM analysis enables to calibrate this new technique to numerical values proposed in the engineering designing codes. Hence, a more accurate probabilistic method may use some experimental-based admissible values included in the existing design of legal provisions.

1. Introduction

Structural safety of civil engineering structures [1] concerning possible fire accidents is a very important practical problem, especially in the area of steel structures, where the demands and expectations are unusually high [2]. Hence, the reliability of engineering structures under fire conditions [3,4] remains always a very challenging and practically important area of knowledge [5,6]. This is due to a partial lack of experimental data, the complexity of the numerical simulation, even while using the Finite Element Method (FEM) or the Finite Volume Method (FVM) as well as the difficulty in stochastic simulations. The full-scale experiments necessary to build and calibrate efficient numerical models [7,8] need to include both realistic fire scenarios, temperature-dependent paths of deformations, and stresses as well as detailed information concerning mechanical and physical characteristics of the structural materials. A fundamental minority in the experimental methods is that individual large-scale structures may be subjected to the fire only once; some structural elements (such as the cold-formed, for instance [9]) or their connections (cf. [10]) may be tested with a few times, but their qualitative and quantitative results have limited application for the framed-structures [11,12], buildings, bridges [13], or steel large-surface halls [14]. It has been widely documented that Computational Fluid Dynamics (CFD) is one of the most powerful numerical simulation tools in fire propagation prediction [15], which is commonly used with some FEM systems to predict the failure [16], and where fire–thermomechanical interface may play a very important role [17].

There is no doubt that deterministic numerical simulations of fire accidents may be very useful in engineering practice, taking into account the dramatic cost of such full-scale experiments, and may have repeatable character; nevertheless, these models are very sensitive to mechanical and thermal boundary conditions, temperature-dependent material parameters, and also the details in their algorithms [18,19]. It can be realized by some hybrid computer systems with an option of the gas–solid interaction or just by the coupled solid mechanics numerical analyses (e.g., thermo-mechanical) available in many well-known FEM commercial systems such as ABAQUS [3,14,19] or ANSYS. Such a simulation could serve as an efficient prediction of the structural failure time to obtain more specific information concerning fire resistance and evacuation time (decisive especially for high-rise buildings safety or large bridges). Let us underline that all these methods and case studies have purely traditional deterministic character and do not allow directly for any structural reliability assessment.

Stochastic models in fire safety analyses are traditionally related to the models of fire spread in woodlands [20], numerical simulation related to this issue [21], as well as to fire outbreaks following earthquake disasters [22]. The Monte-Carlo simulation approach was traditionally the first numerical technique to analyze structural response under a fire [23], and to deliver risk analysis for steel beams [24]. It is well-known that the widely accessible and relatively easily programmable Monte-Carlo simulation needs enormously huge computer time and power consumption [25]; a semi-analytical technique availability depends upon the initial choice of the input uncertainty type [3,26,27], whereas various expansion techniques (Karhunen-Loeve and polynomial chaos [28] or Taylor [26]) may exhibit limited applicability in terms of the larger initial uncertainty level or multivariable (and state-dependent) character of the material and physical characteristics. Some Stochastic Finite Element Method (SFEM) studies are available in the literature [29,30,31], but their connection with fire simulation is rather scarce [3], so no well-documented experiments and related conclusions can be found. The most popular engineering tool in this area was the Second Order Second Moment (SOSM) [32,33], which has been further generalized to the higher order approach engaging the Least Squares Method (LSM) [34,35] recovery of polynomial bases [36] relating the desired structural response with a given input of random parameter(s); it is in fact similar to the response surface methodology [37,38]. The main goal of such an approach would be a final calculation of the reliability index, which in civil engineering designing codes is still based upon the First Order Reliability Method (FORM) [1,39]. Despite the stochastic computer method chosen, the main difficulty would be the collection of basic statistical parameters and corresponding probability distributions for material/physical characteristics of structural materials subjected to high temperatures. It is known that some alternative stochastic approach could be based on a probabilistic distance or probabilistic divergence, but their application is still not quite straightforward in reliability assessment. One of the alternatives in this area could be the so-called Bhattacharyya divergence [14,27,40], but many other models can be useful including Shannon, Renyi, or Tsallis entropies [41,42] and probabilistic distances [43] including Hellinger theory [44], Jeffreys model [45], Kullback–Leibler theory [46], and also Jensen–Shannon entropy [47]. Let us note that certain engineering uncertainty analyses have been delivered in the context of various entropies in the literature [48,49,50,51], but they are a little bit distant from engineering reliability index determination.

The main aim of this paper is to present some stochastic numerical analysis schemes of the fire scenario and to apply them to analyze the reliability of some popular steel structures of the hot-rolled I beam being a part of the structural roof. A very important aspect of this model is its fully coupled character, capturing of temperature variations of all material parameters, and also the usage of 3D finite elements, which enables decisively higher numerical accuracy than the Euler–Bernoulli, Timoshenko, or shell elements applicable in engineering practice; some sequential coupling with ABAQUS has been demonstrated by the authors in [3] before. The stochastic scheme is based upon a triple stochastic methodology–with (i) Monte-Carlo simulation, (ii) semi-analytical approach as well as (iii) iterative generalized higher (the 10th) order stochastic perturbation technique. The entire computational implementation has been carried out with the use of the FEM system ABAQUS^® as well as the computer algebra system MAPLE 2019.2^® and has a general character, independent of the engineering structure type. Determination of the first four probabilistic moments and coefficients available in this algorithm allows for a calculation of the reliability index based on the First Order Reliability Method (FORM), and this index is presented as a function of the fire duration time (and its mean temperature). Additionally, a concept of the probabilistic divergence (relative entropy) usage to approximate structural reliability has been presented and discussed here by a contrast of this entropy to the FORM index [3]. This concept has been successfully employed before to study the reliability of statically uploaded linear elastic steel truss [27], and also in the stochastic dynamic response of some steel halls [14]. The most creative work is an extension of the entropy-based approach from elastic problems towards fully coupled thermo-elasticity Stochastic Finite Element Method structural analysis. Contrary to the previous studies, input uncertainty in fire gas temperatures induces multiple random variability in mechanical and thermal characteristics of the given structure. A methodology proposed here may serve for further fire (and not only) structural safety analysis of the steel structures and the very important aspect of this study is application of the 3D finite elements for detection of structural behavior of steel thin-walled elements. The authors have invented the new reliability index, which reflects probabilistic divergence in-between admissible and extreme deformations depending both on higher temperatures. It has been shown that probabilistic entropy may be efficiently used in engineering analyses not only in the context of the maximum entropy principle but also as a direct function of the input uncertainty and may contain key information concerning structural reliability. The essential innovative aspect of this work is to apply the relative entropy-based reliability index in fully coupled thermo-elastic FEM analysis for simulation of fire accidents in some popular civil engineering structures.

2. Physical Model and Its Implementation

Let us consider a transient thermo-mechanical boundary initial value problem defined on isotropic and homogeneous domain Ω shown schematically below in Figure 1.

Its mechanical part is driven by the following incremental static equilibrium equations [52]:

$$\Delta {\sigma }_{kl,l}+\rho \Delta f_k=0;x\in \mathrm{\Omega }$$

(1)

$$\Delta {\sigma }_{kl}=C_{klmn}\left(T\right)\Delta {\epsilon }_{mn}+{\delta }_{kl}\alpha \left(T\right)\Delta TE\left(T\right);x\in \mathrm{\Omega }$$

(2)

$$\Delta {\epsilon }_{mn}=\frac{1}{2}[\Delta u_{k,l}+\Delta u_{l,k}+u_{i,k}\Delta u_{i,l}+\Delta u_{i,k}{u}_{i,l}+\Delta u_{i,k}\Delta u_{i,l}];x\in \mathrm{\Omega }$$

(3)

with the following essential and natural boundary conditions:

$$\Delta {\sigma }_{kl}{n}_l=\Delta {\widehat{t}}_k;x\in \partial {\mathrm{\Omega }}_{\sigma },k=1,2,3$$

(4)

$$\Delta u_k=\Delta {\widehat{u}}_k;x\in \partial {\mathrm{\Omega }}_u,k=1,2,3$$

(5)

This problem is solved for the displacement vector $u_k\left(x\right)$, the strain tensor ${\epsilon }_{kl}\left(x\right)$ and the stress tensor ${\sigma }_{kl}\left(x\right)$, where the stress tensor increments $\Delta {\sigma }_{kl}\left(x\right)$, $\Delta {\tilde{\sigma }}_{kl}\left(x\right)$ and denote the first and the second Piola–Kirchhoff tensors

$$\Delta {\sigma }_{kl}=\Delta F_{km}\Delta {\tilde{\sigma }}_{ml}+F_{km}\Delta {\tilde{\sigma }}_{ml}+\Delta F_{km}{\tilde{\sigma }}_{ml};x\in \mathrm{\Omega }$$

(6)

with

$$\Delta F_{km}=\Delta u_{k,m};x\in \mathrm{\Omega }.$$

(7)

All static state variables, i.e., displacements, strains and stresses are temperature-dependent, but this dependence is omitted for a brevity of presentation in all equilibrium equations. Simultaneously, a transient heat flow problem for the temperature field T = T(x,τ) is solved from the following differential equation [53]:

$$\rho \left(T\right)c\left(T\right)\dot{T}-{\left({\lambda }_{ij}\left(T\right)T_{,j}\right)}_{,i}-g=0;x_i\in \mathrm{\Omega };\tau \in [0,\infty ),$$

(8)

where c(T) is the temperature-dependent heat capacity of the region Ω, ρ(T) is the temperature-dependent material density of Ω, λ_ij(T) is temperature-dependent second-order tensor thermal conductivity, and g is the rate of heat generated per unit volume. Traditionally, T and τ denote temperature field values and time, respectively.

This equation should fulfil the boundary conditions of the ∂Ω, which are given as follows:

(1)

temperature (essential) boundary conditions

$$T=\widehat{T};x\in \partial {\mathrm{\Omega }}_T,$$

(9)

and for ∂Ω_q part of the total ∂Ω:

(2)

heat flux (natural) boundary conditions

$$\frac{\partial T}{\partial n}=\widehat{q};x\in \partial {\mathrm{\Omega }}_q,$$

(10)

where $\partial {\mathrm{\Omega }}_T\cup \partial {\mathrm{\Omega }}_q=\partial \mathrm{\Omega }$ and $\partial {\mathrm{\Omega }}_T\cap \partial {\mathrm{\Omega }}_q=\{\varnothing \}$. The initial conditions have been introduced as

$$T^0=T\left(x_i;0\right);x_i\in \mathrm{\Omega },\tau =0.$$

(11)

The following functional defined on $\Delta u_k$ is introduced to obtain a numerical solution to the deformation problem:

$$J\left(\Delta u_k\right)={\displaystyle \underset{\mathrm{\Omega }}{\int }\left(\frac{1}{2}{C}_{klmn}\left(T\right)\Delta {\epsilon }_{kl}\Delta {\epsilon }_{mn}+\frac{1}{2}{\tilde{\sigma }}_{kl}\Delta u_{i,k}\Delta u_{i,l}-\rho \Delta f_k\Delta u_k\right)d\mathrm{\Omega }}-{\displaystyle \underset{\partial \mathrm{\Omega }}{\int }\Delta {\widehat{t}}_k\Delta u_{k}d\left(\partial \mathrm{\Omega }\right)}$$

(12)

whose solution is determined from the minimization of the incremental version of the potential energy stationarity principle

$$\delta J\left(\Delta u_k\right)=\frac{\partial J}{\partial \Delta u_k}\delta \left(\Delta u_k\right)$$

(13)

Analogously, one considers some continuous temperature variations $\delta T\left(x_i\right)$ defined in the interior of the region Ω and vanishing on $\partial {\mathrm{\Omega }}_T$. A variational formulation may be proposed here as

$$\begin{array}{c}{\displaystyle \underset{\mathrm{\Omega }}{\int }\left(\rho \left(T\right)c\left(T\right)\dot{T}\delta T+{\lambda }_{ij}^{}\left(T\right)T_{,j}\delta T_{,i}-g\delta T\right)}d\mathrm{\Omega }-{\displaystyle \underset{\partial {\mathrm{\Omega }}_q}{\int }\widehat{q}\delta Td\left(\partial \mathrm{\Omega }\right)}=0;\\ x_i\in \mathrm{\Omega };\tau \in [0,\infty ).\end{array}$$

(14)

Let us recall the classical Finite Element Method formulation, where the displacements increments $\Delta u_k\left(x\right)$ being a continuous and differentiable function over the region Ω consisting of the geometrically continuous subsets (finite elements) ${\mathrm{\Omega }}_e$, where e = 1,…,E gives a complete representation of the set Ω. Let us consider the following approximation of the displacement increments [52,54]:

$$\Delta u_k\left(x\right)={\displaystyle \sum _{\zeta =1}^{N_e}{\varphi }_{\zeta k}\left(x\right)}\Delta q_{\zeta }^{(N)},$$

(15)

where ${\varphi }_{\zeta k}\left(x\right)$ are the shape functions in the node k, $\Delta q_{\zeta }^{(N)}$ represent the nodal degrees of freedom vector, while N_e is the total number of those degrees of freedom in the considered node. Starting from the proposed approximation it is possible to express the gradients of the displacement vector as well as the strain tensor components as

$$\Delta u_{k,l}\left(x\right)={\varphi }_{k,l}^{\zeta }\left(x\right)\Delta q_{\zeta }^{(N)},$$

(16)

$$\Delta {\overline{\epsilon }}_{kl}\left(x\right)=[{\overline{B}}_{kl}^{(1)\zeta }+{\overline{B}}_{kl}^{(2)\zeta }]\Delta q_{\zeta }^{(N)}={\overline{B}}_{kl}^{\zeta }\Delta q_{\zeta }^{(N)},$$

(17)

$$\Delta {\overline{\overline{\epsilon }}}_{kl}\left(x\right)={\overline{\overline{B}}}_{kl}^{\zeta \xi }\Delta q_{\zeta }^{(N)}\Delta q_{\xi }^{(N)}$$

(18)

and finally

$$\Delta {\epsilon }_{kl}\left(x\right)=\Delta {\overline{\epsilon }}_{kl}\left(x\right)+\Delta {\overline{\overline{\epsilon }}}_{kl}\left(x\right).$$

(19)

The following notation has been applied in the above equations:

$${\overline{B}}_{kl}^{(1)\zeta }\left(x\right)={\varphi }_{k,l}^{\zeta }\left(x\right),{\overline{B}}_{kl}^{(2)\zeta }\left(x\right)={\varphi }_{i,k}^{\zeta }\left(x\right){\varphi }_{i,l}^{\xi }\left(x\right)q_{\xi }^{(N)},{\overline{\overline{B}}}_{kl}^{\zeta \xi }\left(x\right)=\frac{1}{2}{\varphi }_{i,k}^{\zeta }\left(x\right){\varphi }_{i,l}^{\xi }\left(x\right).$$

(20)

Now, the following elemental stiffness matrices are introduced

$$k_{\zeta \xi }^{(\sigma )e}={\displaystyle \underset{{\mathrm{\Omega }}_e}{\int }{\tilde{\sigma }}_{kl}{\varphi }_{i,k}^{\zeta }\left(x\right){\varphi }_{i,l}^{\xi }\left(x\right)d\mathrm{\Omega }},$$

(21)

$$k_{\zeta \xi }^{(con)e}={\displaystyle \underset{{\mathrm{\Omega }}_e}{\int }\frac{1}{2}{C}_{klmn}\left(T\right){\overline{B}}_{kl}^{(1)\zeta }{\overline{B}}_{mn}^{(1)\xi }d\mathrm{\Omega }},$$

(22)

$$k_{\zeta \xi }^{(u)e}={\displaystyle \underset{{\mathrm{\Omega }}_e}{\int }\frac{1}{2}{C}_{klmn}\left(T\right)\left({\overline{B}}_{kl}^{(1)\zeta }{\overline{B}}_{mn}^{(2)\xi }+{\overline{B}}_{kl}^{(2)\zeta }{\overline{B}}_{mn}^{(1)\xi }+{\overline{B}}_{kl}^{(2)\zeta }{\overline{B}}_{mn}^{(2)\xi }\right)d\mathrm{\Omega }},$$

(23)

while the second and the third order stiffnesses are equal, respectively

$$k_{\zeta \xi }^{(2)e}={\displaystyle \underset{{\mathrm{\Omega }}_e}{\int }\frac{3}{2}{C}_{klmn}\left(T\right)\left({\overline{B}}_{kl}^{\zeta }\Delta u_{\zeta }^{(N)}{\overline{\overline{B}}}_{mn}^{\mu \nu }\Delta q_{\mu }^{(N)}\Delta q_{\nu }^{(N)}+{\overline{\overline{B}}}_{kl}^{\zeta \xi }\Delta q_{\zeta }^{(N)}\Delta q_{\xi }^{(N)}{\overline{B}}_{mn}^{\mu }\Delta q_{\mu }^{(N)}\right)d\mathrm{\Omega }},$$

(24)

$$k_{\zeta \xi }^{(3)e}={\displaystyle \underset{{\mathrm{\Omega }}_e}{\int }2C_{klmn}\left(T\right)\left({\overline{\overline{B}}}_{kl}^{\zeta \xi }\Delta q_{\zeta }^{(N)}\Delta q_{\xi }^{(N)}{\overline{\overline{B}}}_{mn}^{\mu \nu }\Delta q_{\mu }^{(N)}\Delta q_{\nu }^{(N)}\right)d\mathrm{\Omega }}.$$

(25)

Introducing $k_{\zeta \xi }^{(i)}$ for i = 1,2,3 into the functional $J\left(\Delta u_k\right)$ in Equation (8), applying a transformation from the local to global coordinates system, one may obtain from the stationarity condition that

$$K_{\alpha \beta }^{(1)}\left(T\right)\Delta q_{\beta }+K_{\alpha \beta \gamma }^{(2)}\left(T\right)\Delta q_{\beta }\Delta q_{\gamma }+K_{\alpha \beta \gamma \delta }^{(3)}\left(T\right)\Delta q_{\beta }\Delta q_{\gamma }\Delta q_{\delta }=\Delta Q_{\alpha }$$

(26)

Fulfilled for any configuration of the region Ω, where ${\widehat{K}}_{\alpha \beta }$, $q_{\beta }$, and $Q_{\alpha }$ are stiffness matrix, displacement vector, and nodal loads vector, respectively. The same discretization serves for the discretization of the temperature field by the nodal temperatures vector ${\theta }_{\alpha }$ as [52,54]

$$T(x_i)=H_{\alpha }\left(x_i\right){\theta }_{\alpha };i=1,2;\mathsf{\alpha }=1,2,\dots ,\mathrm{N},$$

(27)

The temperature gradients can be rewritten in the form

$$T_{,i}=H_{\alpha ,i}{\theta }_{\alpha },\mathrm{i}=1,2.$$

(28)

We introduce the capacity matrix $C_{\alpha \beta }\left(T\right)$, the heat conductivity matrix ${\widehat{K}}_{\alpha \beta }\left(T\right)$ and the vector $P_{\alpha }$ as

$$C_{\alpha \beta }\left(T\right)={\displaystyle \underset{\mathrm{\Omega }}{\int }\rho \left(T\right)c\left(T\right)H_{\alpha }{H}_{\beta }d\mathrm{\Omega }},\hspace{1em}\hspace{1em}{\widehat{K}}_{\alpha \beta }\left(T\right)={\displaystyle \underset{\mathrm{\Omega }}{\int }{\lambda }_{ij}^{}\left(T\right)H_{\alpha ,i}{H}_{\beta ,j}d\mathrm{\Omega }}$$

(29)

and also the R.H.S. vector in the following way:

$$P_{\alpha }={\displaystyle \underset{\mathrm{\Omega }}{\int }g{H}_{\alpha }d\mathrm{\Omega }}+{\displaystyle \underset{\partial \mathrm{\Omega }}{\int }\widehat{q}}{H}_{\alpha }d\mathrm{\Omega }.$$

(30)

Next, let us introduce these matrixes into the variational formulation (14) to obtain the following algebraic equations system:

$$C_{\alpha \beta }\left(T\right){\dot{\theta }}_{\beta }+{\widehat{K}}_{\alpha \beta }\left(T\right){\theta }_{\beta }=P_{\alpha }$$

(31)

Equations (26) and (31) are finally solved simultaneously by the system ABAQUS to obtain time fluctuations of the deformations and stresses into the given boundary thermo-elasticity problem. This coupled problem is solved in a homogeneous continuous domain with no initial stresses and strains. A steel material occupying this domain shows temperature-dependent material and physical characteristics presented in Figure 2 and Figure 3 below. It is reported here, after experimental works included in Eurocode 3 [2], that Young modulus as well as heat conductivity decrease while increasing temperature induced into the steel volume, quite opposite to thermal elongation, which increases. Heat capacity shows very complex behavior in this situation, especially in terms of the singularity reached close to 700 °C. This singularity results in remarkable difficulty during a solution of Equation (31), which is practically omitted by listing heat capacity at each 50 °C and numerical interpolation in-between these discrete data.

3. Probabilistic Aspects and Relative Entropy

Let us introduce the random variable b and its probability density function as $p_b\left(x\right)$. Then, the first two probabilistic moments of this variable are defined as [55]

$$E\left[b\right]\equiv b^0={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}b{p}_b\left(x\right)dx}$$

(32)

where $b^0$ means the average value of $b$ itself and

$$Var\left[b\right]={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\left[b-E\left(b\right)\right]}^2{p}_b\left(x\right)dx}$$

(33)

Higher probabilistic moments and related coefficients may be defined according to the classical definitions from the probability theory. The basic idea of the stochastic perturbation approach employed here is to expand all input random variables and all the resulting state functions of the given boundary initial problem via the Taylor series about their spatial expectations using the perturbation parameter $\epsilon$. The random function $f\left(b\right)$ with respect to its parameter $b$ about its mean value is given as follows

$$f=f^0+{\displaystyle \sum _{n=10}^1\frac{{\epsilon }^n}{n\mathit{!}}}\frac{{\partial }^{n}f}{\partial b^n}\Delta b^n,\mathrm{where}\epsilon \Delta b=\epsilon \left(b-b^0\right)$$

(34)

is the first variation of a variable $b$ about its expected value and the symbol ${(.)}^0$ represents the value of a function calculated for its mean. Let us analyze further the expected values of displacement state function $u\left(b\right)$ defined by its expansion via the Taylor series as follows

$$E\left[u\left(b\right)\right]={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}u\left(b\right)p_b\left(x\right)dx}={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\left(u^0\left(b^0\right)+{\displaystyle \sum _{n=10}^1\frac{{\epsilon }^n}{n\mathit{!}}}\frac{{\partial }^{n}u\left(b\right)}{\partial b^n}\Delta b^n\right)}{p}_b\left(x\right)dx$$

(35)

This expansion is valid only if the state function is analytic in $\epsilon$ and the series converges and, therefore, any criteria of convergence should include the magnitude of the perturbation parameter. A perturbation parameter here as usually in practical engineering computations is equal 1. It yields for the input random variable with symmetric probability density function in the tenth order approach

$$E\left[u\left(b\right)\right]=u^0\left(b^0\right)+\frac{1}{2\mathit{!}}\frac{{\partial }^{2}u\left(b\right)}{\partial b^2}{\mu }_2\left(b\right)+\cdots +\frac{1}{10\mathit{!}}\frac{{\partial }^{10}u\left(b\right)}{\partial b^{10}}{\mu }_{10}\left(b\right)$$

(36)

${\mu }_n\left(b\right)$ denotes the nth order central probabilistic moment of variable b. All the terms with odd orders are equal to 0 for the symmetric random variable and the orders higher than the 10th have been simply neglected. Similar considerations lead to the expressions, such as the variance, for instance

$$Var\left[u\left(b\right)\right]={\mu }_2\left(u\left(b\right)\right)={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\left(u\left(b\right)-E\left[u\left(b\right)\right]\right)}^{2}p\left(b\right)db}$$

(37)

Quite similarly, it is possible to derive higher-order central probabilistic moments from their definitions

$${\mu }_3\left(u\left(b\right)\right)={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\left(u\left(b\right)-E\left[u\left(b\right)\right]\right)}^{3}p\left(b\right)db},\hspace{1em}\hspace{1em}{\mu }_4\left(u\left(b\right)\right)={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\left(u\left(b\right)-E\left[u\left(b\right)\right]\right)}^{4}p\left(b\right)db}$$

(38)

Let us mention that it is necessary to multiply each of these equations by the relevant order probabilistic moments of the input random variable to obtain the algebraic form convenient for any symbolic computations. Based on the classical definition of the variance we can calculate the coefficient of variation, skewness, and kurtosis as follows

$$\alpha \left(u\left(b\right)\right)=\frac{\sqrt{{\mu }_2\left(u\left(b\right)\right)}}{E\left[u\left(b\right)\right]},S\left(u\left(b\right)\right)=\frac{{\mu }_3\left(u\left(b\right)\right)}{{\left[\sqrt{{\mu }_2\left(u\left(b\right)\right)}\right]}^3},\kappa \left(u\left(b\right)\right)=\frac{{\mu }_4\left(u\left(b\right)\right)}{{\left[\sqrt{{\mu }_2\left(u\left(b\right)\right)}\right]}^4}-3$$

(39)

It should be mentioned that at this stage the proposed procedure is still independent of a choice of the initial probability distribution function, however, a satisfactory probabilistic convergence of the final results may demand various lengths of the expansions for random variables with different distributions. A common implementation of the Monte-Carlo simulation, semi-analytical approach, as well as higher order stochastic perturbation technique in the system MAPLE was possible thanks to the Least Squares Method polynomial recovery of the structural displacements with respect to the input uncertain parameter b.

Finally, the reliability index is to be determined to measure the structural safety of the given case study of the beam under fire conditions. The engineering codes (such as Eurocode 0) advise applying the FORM theory, where one calculates

$${\beta }_{FORM}=\frac{\mathrm{E}\left[R\right]-\mathrm{E}\left[E\right]}{\sqrt{Var\left(R\right)+Var\left(E\right)}};$$

(40)

E[R] denotes here the maximum acceptable displacement of a beam according to Eurocode 3, and E[E] is the expected value of extreme displacement calculated based upon the iterative generalized stochastic perturbation technique. It is well known that the FORM approach has some limitations, so that higher order theories have been developed such as the Second and even Third Order Reliability Methods (SORM and TORM), especially because of the linear character of the limit function in the FORM technique. This seems to be unjustified in many practical applications and undoubtedly fire simulation is one of them. Further, the main mathematical methodological difficulty is that both distributions, of R and of E, are assumed as Gaussian here, which may be far from the experimental evidence and even misleading in the view of many theoretical studies and computer simulations.

This was a reason to seek for another concept to assess structural reliability in the presence of fire and to define the reliability index using probabilistic distance in-between two distributions of structural resistance and of structural extreme effort. One of the available mathematical models is the Bhattacharyya theory [40], which proposes for two different PDFs, namely p_R(x) and p_E(x), the following probabilistic distance measure:

$${\beta }^{\prime }={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\sqrt{p_R(x)p_E(x)}dx},$$

(41)

where associated with structural resistance and probability function related to structural effort, respectively. Such formula enables the application of two different probability distributions of practically any nature and the given set of their parameters; further application of the analytical integration using some computer algebra system may lead to the desired numerical value. However, this formula can be simplified in the case of two Gaussian distributions with the given expectations and standard deviations (μ(R), μ(E), and also σ(R), σ(E), correspondingly). There holds

$${\beta }^{\prime }=\frac{1}{4}\frac{\left(E\left[R\right]-E\left[E\right]\right)}{{\sigma }^2\left(R\right)+{\sigma }^2\left(E\right)}+\frac{1}{2}\ln \left(\frac{{\sigma }^2\left(R\right)+{\sigma }^2\left(E\right)}{2\sigma \left(R\right)\sigma \left(E\right)}\right)$$

(42)

Another problem is the upscaling of the variability interval of reliability indices obtained as the result of Equation (41) (or (42)) to contrast them with the values resulting from the FORM approach. The main idea would be to apply this mathematical apparatus with its theoretical importance, but to retain the existing engineering codes recommended minimum values of the reliability indices. Without such an upscaling direct determination of the reliability index, it would never show any extremes separating safe and unsafe design and/or exploitation of the given structure. Some preliminary computational experiments enable to propose the following rescaling of such a reliability index to numerical values comparable with the existing codes:

$$\beta =\frac{1}{2}\sqrt{\frac{1}{4}\frac{\left(E\left[R\right]-E\left[E\right]\right)}{{\sigma }^2\left(R\right)+{\sigma }^2\left(E\right)}+\frac{1}{2}\ln \left(\frac{{\sigma }^2\left(R\right)+{\sigma }^2\left(E\right)}{2\sigma \left(R\right)\sigma \left(E\right)}\right)}$$

(43)

Both square root function as well as the multiplier preceding this function has been taken to make the first component, free from logarithmic function as close as possible to Equation (40); further numerical analysis would check an efficiency of such a modified reliability index. Computational implementation of the aforementioned approach has hybrid character and uses both the FEM system ABAQUS Standard, and also computer algebra package MAPLE. ABAQUS Standard fully coupled thermoelastic analysis is run multiple times to have the series of temperatures, stresses, and displacements obtained for the few different input fire temperatures. Then, these discrete results are in vector formats transferred to the MAPLE system, where the Least Squares Method procedures enable to recover polynomial approximations of these state functions with respect to fire temperature. Finally, this fire temperature is randomized according to the Gaussian distribution with the given mean value and given variability interval of the coefficient of variation. Polynomial approximations of temperatures, stresses, and displacements together with a triple implementation of probabilistic analysis return the basis probabilistic characteristics (i.e., expectations and standard deviations) of the state functions. Finally, two different reliability indices presented in Equations (40) and (43) are determined as the functions of fire exposure time.

4. Computational Illustration

The numerical example consists of the 6m both ends fixed steel I-profile beam with linear load applied to top flange equal 2.0 kN/m (Figure 4), whereas its cross section and the FEM discretization have been shown in Figure 5. The mesh consisting of 34,800 brick finite elements has been used; they are denoted in the system ABAQUS as C3D8T: the 8-node thermally coupled brick, tri-linear displacement, and temperature. The edge length of each brick element is less than 10 mm and it guarantees quite smooth approximation of stress distribution throughout the web. The thermal load has been adopted from the Standard ISO fire curve (Figure 6). According to Eurocode fire fumes after time t_f = 15 min are about 740 ºC. This thermal load has been applied to each side of the cross-section except the top flange where the concrete slab lays according to the given fire scenario (Figure 6). Surface radiation and surface film conditions as thermal loads have been applied in this model (Figure 7). All material parameters are fully temperature-dependent and the fully coupled transient temperature-displacement incremental analysis has been used with the full Newton solution technique. The total time period of numerical simulation has been set as 900 s. The initial increment size is 0.01 step time, in this case, the minimum increment size is 0.0001 of the step time, and the maximum increment has been fixed as a single second.

Figure 8 shows a comparison of this beam deflection under the given fire conditions (left diagram) and without a fire (right diagram). It is noticeable that fire heating increases the extreme deflection at half of the beam span more than four times (13.0 cm) than its admissible value (3.0 cm); maximum deflection in this case of no fire is less than 2.0 mm. One concludes that extreme deformation in the beam under fire is two orders larger than for the beam with no temperature attached (analytical approximation following basic equations of strength of materials appears to be quite efficient). The second major difference in-between these two models is that a beam subjected to fire shows remarkable deformations close to the supported area, whereas the beam without fire loadings has these deformations negligibly small. Next, Figure 9 reports normal stress marked as S33 (left diagram) and the reduced von Mises stress (right diagram). Von Mises stresses are almost 4 times larger than the longitudinal stresses in this specific case and both stress fields exhibit quite different distributions along the beam and throughout its cross-section. Figure 10 shows for a completeness temperature distribution at the steady-state (at the end of the beam heating by a fire) and one compares with an engineering intuition that the thinner web accumulates definitely more heat than the relatively thicker rest of the beam.

Further results concern probabilistic characteristics of the structural response, which include expectations, coefficients of variation, skewness, and kurtosis of displacements shown in Figure 11, Figure 12, Figure 13 and Figure 14, correspondingly. They have been computed using three different probabilistic methods, namely the perturbation method (abbreviated as a PM here), the semi-analytical method (SAM), and also using the Monte-Carlo simulation (MCS); they have been shown all as the functions of the input coefficient of variation ranging from 0.0 until 0.25. They concern extreme vertical displacements obtained at the half of the beam structure for the needs of further reliability assessment according to the Serviceability Limit State (SLS). First of all, it is seen that the first two moments are stable until t = 500 s, and then they start to diverge; initially, both moments equal almost 0. Generally, three different probabilistic numerical methods coincide in the case of the first two probabilistic characteristics with each other until α = 0.20. Higher order statistical characteristics are more dispersed—they exhibit both numerical values very close to 0, which start to diverge at about t = 500 s. This can be interpreted that until this time moment of the fire exposure, the resulting displacements can be approximated as Gaussian, which means that both the FORM expression for the reliability index as well as its relative entropy counterpart proposed by Equation (43) are justified very well. Finally, one can notice that all three probabilistic methods coincide very well until α(t_f) = 0.10.

The next set of probabilistic characteristics computed for the reduced von Mises stresses is given similarly in Figure 15, Figure 16, Figure 17 and Figure 18. They have been collected here as decisive in the calculation of the reliability index in the Ultimate Limit State (ULS). The expected values increase monotonously from 0 until the same moment t = 500 s; then, they exhibit large fluctuations, although all three probabilistic methods return the same numerical values. The reduced von Mises stresses start to highly depend upon the input coefficient of variation for the fire exposure time t > 500 s, which is an observation quite unusual for elastic problems with any uncertainty. A little bit different conclusions can be drawn from the results contained in Figure 16—the output CoVs keeps very close to 0 during the entire fire heating process showing some numerical discrepancies almost at the beginning of this process (t = 50 s) as well as its end (t = 800 s). Higher order statistics are rather distant from 0 (see Figure 17 and Figure 18), so they cannot be efficiently modeled as Gaussian and need larger numerical effort. A coincidence of Monte-Carlo simulation, the semi-analytical method as well as the stochastic perturbation technique is worse and can be assumed α(t_f) = 0.05. This is a quite expected result because probabilistic characteristics of the stresses are calculated based on probabilistic moments of displacements (since displacement-version of the FEM is used) and of probabilistic characteristics of the constitutive tensor (additionally depending on the nodal temperatures).

Higher order statistics of the reduced stresses are definitely distant from 0 during fire exposure time, so that these stresses cannot be approximated with the Gaussian distribution. Therefore, the FORM index given in Equation (40) includes a remarkable modelling error, whereas relative entropy should be calculated here using the general analytical formula (42). Further numerical results concern reliability analysis of the given structure, so that Figure 19 contains the key results in this study and presents a comparison of the reliability indices determined using the FORM approach (left graph) with these calculated thanks to the Bhattacharyya relative entropy. They are both plotted as the functions of the fire exposure time and are related to the Serviceability Limit State. This is because previous results document well approximation of the distribution of the displacement by the Gaussian PDF; this is not the case with the reduced stresses, so their reliability has been discussed in Figure 20 using the FORM index only. As one could expect, these time fluctuations of reliability indices have typical exponential decay in the first stage of fire heating, but then, after reaching critical and 0 values, exhibit further oscillatory character. This part is interesting from the numerical point of view and has no practical importance and the structure lost its reliability. Such a wavy character is totally absent for civil engineering structures subjected to, i.e., static loads and some stochastic ageing process, when reliability decreases monotonously.

The very important conclusion following the comparison of the two graphs in Figure 19 is that the FORM index and relative entropy variations show the same character and extreme values location. The most important conclusion is that they cross the admissible value at almost the same time as a fire accident, which means that structural safety analysis results in both cases in the same evacuation time. Further, it is seen that the FORM analysis results in a few times higher reliability index at the very beginning of the fire accident. This can be meant as some overestimation of the realistic reliability, but this difference is meaningless for structural design. A very interesting result is that the FORM index is more sensitive to the input uncertainty level, whereas the relative entropy shows almost no such sensitivity.

The reliability index in the ULS analysis shows less waviness above the reliability limit. It does not affect the conclusion and structural safety—it concerns computational aspects rather. A contrast of Figure 19 and Figure 20 results in a conclusion that the ULS is decisive for the overall safety of this element—it seems that the given beam faster falls into the plastic regime than approaches the admissible deformations at half of its span. This conclusion is a little bit out of the probabilistic analysis, nevertheless, the stochastic reliability study confirms an engineering observation. This notice is supported mainly by the reduced von Mises stress time fluctuations. Evaluation of the beam safety from the normal longitudinal stresses may lead to an improper conclusion that the ULS and the SLS exhibit almost the same failure time.

A more detailed comparison of the classical FORM approach with the proposed new one based upon the relative entropy has been provided in

Table 1. Reliability index for the vertical displacements as a function of fire duration values comparison for the FORM theory and the relative entropy approach.

	βFORM(uz(tf))					β(uz(tf))
tf [sec.]	α(tf)					α(tf)
	0.05	0.10	0.15	0.20	0.25	0.05	0.10	0.15	0.20	0.25
0	1378.65	1379.21	1367.60	1369.04	1371.28	1378.93	1378.87	1378.81	1378.75	1378.65
60	312.18	508.92	397.98	359.21	340.79	312.26	312.29	312.34	312.32	312.18
120	−97.00	130.42	10.54	−34.75	−57.39	−96.87	−96.84	−96.81	−96.85	−97.00
180	−153.26	23.50	−63.06	−98.32	−116.80	−152.95	−152.98	−153.04	−153.12	−153.26
240	−65.15	36.38	−6.63	−26.80	−38.25	−64.67	−64.79	−64.91	−65.03	−65.15
300	38.66	74.53	66.69	59.70	54.77	39.26	39.07	38.90	38.76	38.66
360	93.80	89.41	100.40	101.46	100.87	94.41	94.19	94.00	93.87	93.80
420	84.44	67.40	80.84	84.37	85.56	84.97	84.74	84.58	84.48	84.44
480	27.98	18.82	23.52	25.38	26.31	28.35	28.16	28.05	28.00	27.98
540	−40.04	−32.90	−40.25	−41.77	−42.08	−39.87	−40.00	−40.05	−40.04	−40.04
600	−80.72	−61.70	−76.53	−80.50	−81.90	−80.75	−80.81	−80.80	−80.75	−80.72
660	−66.04	−49.14	−61.33	−64.89	−66.32	−66.24	−66.24	−66.19	−66.10	−66.04
720	7.35	5.74	7.06	7.25	7.16	7.06	7.09	7.15	7.26	7.35
780	104.99	77.01	96.39	101.92	103.94	104.74	104.71	104.78	104.88	104.99
840	142.07	102.81	129.23	137.05	140.10	142.02	141.82	141.89	141.98	142.07
900	−28.24	−23.07	−27.43	−28.62	−29.03	−27.94	−28.46	−28.37	−28.29	−28.24

below (β_FORM(u_z(t_f))) and relative entropy approach (β(u_z(t_f))). Both reliability indices have been compared with each other throughout the entire fire accident simulation time for the few different values of input coefficient of variation. There is no doubt that both methods return almost the same values (positive and negative also) for any input statistical scattering. Furthermore, a failure time defined as the moment of fire duration when the reliability index falls down below the admissible values suggested in the designing codes (t = 80 s) in both methods is also almost the same. It seems that the rescaled relative entropy calculated from Equation (43) enables to predict fire safety with the use of the existing engineering codes.

The data presented in both Table 1, and also in Figure 19 and Figure 20 demonstrate that steel beams without any fire protection are extremely sensitive to fire temperatures. A sensitivity of its deformation is of course a few times higher than of the ULS limit function, which is confirmed by the negative reliability indices detected in case of the SLS. It is important that both reliability approaches return the same qualitative results.

Finally, it should be noticed that quite a satisfactory accuracy of both probabilistic simulations and reliability analysis shows that the methodology presented could be applied to other types of steel structures and further applications towards aluminum alloys may be taken into account. Undoubtedly, this approach would be closer to industrial applications when the phase change (from solid to liquid) could be accounted for; this needs brand new implementations in the system ABAQUS. It should be mentioned that the civil engineering designing codes still do not include any statements enabling efficient engineering reliability concerning fire safety (neither for steel nor for concrete or traditional structures). On the other hand, closer interoperability of the FEM system ABAQUS, with the computer algebra packages such as MAPLE should be also achieved.

5. Conclusions

(1)

The major original contribution of this study is a development of new reliability assessment procedure for steel structures subjected to fire conditions. It has been completed with the use of the relative probabilistic entropy model developed by Bhattacharyya and contrasted after some numerical rescaling with the FORM index inherent in the current designing civil engineering codes. The key research finding in this work is that the Bhattacharyya relative entropy sensitivity with respect to the given fire temperature curve (and the additional duration time) is the same as in the case of the FORM reliability index. The innovative scaling procedure enables to use the new mathematical theory to verify the structural safety according to the existing procedures and recommended values of the reliability indices. This is exactly the same scaling as in the case of linear elastic models of steel trusses and spatial structures of the halls [14,20], which means that this procedure may have general character. It should be mentioned that such a new verification procedure may be relatively easily implemented in any hybrid computational FEM and computer algebra system numerical environment.

(2)

A triple probabilistic approach based on the Monte-Carlo simulation, generalized iterative stochastic perturbation technique, and also the semi-analytical approach programmed in conjunction with the polynomial bases of the fire temperature returns very consistent and stable numerical results of expected values, coefficients of variations, as well as skewness and kurtosis—all as the functions of the fire exposure time. Application of such a technique enabled to discover that extreme deformations of the given system have a probability distribution very similar to the Gaussian one, and hence, determination of their probabilistic divergence could proceed with the use of the analytical formula. Probabilistic divergence of the admissible and extreme stresses needs to be proceeded by semi-analytical integration of the general formula. Despite the rather simple character of the structural beam investigated in this work, such a coincidence of three different probabilistic techniques is expected for other structures and their fire scenarios, which deserves further numerical simulations.