Numerical model for determining fire behaviour of structures

A numerical model and computer program for predicting behaviour of structures subjected to fire action are presented in the paper. The nonlinear numerical procedure is conducted in pre-defined time increments. At that, the distribution of temperature is calculated in each increment and, depending on this calculation, material properties and stiffness of the element are corrected, and the static problem is resolved. The efficiency and accuracy of the model and computer program are presented on an example of a simply supported beam.


Autori:
Izvorni znanstveni rad Original scientific paper Wissenschaftlicher Originalbeitrag Numerical model for determining fire behaviour of structures

Introduction
Numerical modelling of behaviour of structures subjected to fire is currently one of topical subjects of vital importance on the global level.The development of simple and efficient numerical models, especially those that have been confirmed through experimental study [1,2], is the basis for better understanding of fire as a stochastic process, and for further development of construction standards with requirements for better and safer structures.
First scientific studies on the effect of fire on man-made structures were conducted already in 1960s.Since then, development of computers has spurred astonishing advances in all fields of engineering, and hence also fast evolution of mathematical/numerical models for the description of fire action on structures.First models were developed in 1970s, mostly for steel structures, and were based on the hybrid model of nonlinear bearing capacity of cross section, as combined with the heat transfer model.These models most often originated from the 2D heat transfer model that provided sufficiently good predictions about development of temperature in structure, if the structure is uniformly heated along the entire span, which was most often not true in real fire situations.
In more recent studies [3,4,5,6], hybrid models have been complemented and improved with various mechanical structural model formulations based on 1D/2D elements for description of structure/section geometry, which was combined with 3D heat transfer models that enable more accurate prediction of temperature distribution within the structure.The definition of constitutive law of behaviour of materials at high temperatures in case of multiaxial stress still remains a considerable problem in experimental research, especially as it does not enable the use of the 3D mechanical model of cross section.The prediction of fire behaviour of structures is more accurate in case of steel structures when compared to concrete structures.This is due to less complicated behaviour of steel at high temperatures, and also to smaller deviation of experimentally defined mechanical properties of steel [7].The prediction of fire behaviour and fire resistance of structures has become more accurate after introduction of various implicit and explicit material creep models defining creep effect at high temperatures [8,9].

Introduction
The numerical model that accurately depicts behaviour of structures under fire must be capable of describing, in addition to nonlinear behaviour of structures subjected to load, the development of temperature within the structure and the change of material properties at high temperatures.In order to conduct such a complete (nonlinear) analysis, it is of high significance to know the cross-sectional geometry, the exact type and position of reinforcement (for reinforcedconcrete structures), load conditions, and the constitutive material law, which is generally nonlinear.The results of this analysis can greatly alter the stress and strain situation in individual structural elements, and enable structural engineer to gain a better insight into the behaviour and possible bearing capacity failure, which ultimately leads towards construction of more durable and economical structures.

Linear elastic model for beam elements
The principal starting point of almost all nonlinear analyses is the linear analysis, i.e. the analysis of materials conforming to Hooke' s law.As this numerical model for behaviour of beamelement structures is very well known, and as it has often been described in the literature [10,11,12], it will only be briefly outlined in this study.Two-node, rectilinear, ideally straight and partly prismatic final elements, with 6 degrees of freedom per node, considered so far in many papers [11,13,14,15,16,17], are used in this paper.Behaviour of a beam element under load can generally be described using the following linear differential equation of equilibrium: Where: -internal force vector -bracing vector -load vector -differential operator In the absence of analytical solutions, the solution to equation ( 1) is normally sought through numerical procedures.One of the most often used and widely recognized procedures is the finite element method.The basic idea behind this method is to replace the system having an indefinite number of degrees of freedom with the system having a finite number of degrees of freedom.To achieve this, the behaviour of a set of points within a system is assumed (programmed) on a single finite element, as related to a specified number of fixed, previously determined points (nodes) on the same element.An approximate solution of the displacement field on a single element is assumed as follows: where H is the basic function matrix, while is the vector of unknown node displacements.
In case of beam-element systems, basic (shaping) functions are most frequently selected for the Hermite polynomial group [11,13].
The following can be deduced from the external and internal force equivalence: that is: i.e. after the left side is multiplied with : or more simply: where: -vector of internal forces at the ends of the finite element -element stiffness matrix -vector of external forces Numerical model for determining fire behaviour of structures Local stiffness matrices and load vector must be transformed to the global coordinate system and, after transformation, the global system equilibrium is established by simple arrangement of stiffnesses and fixity forces in the corresponding nodes of the finite element mesh.
These are generally known equations where K and F are matrices of stiffness and load, respectively, while u is the global displacement vector.Before solving the above equation system, it is necessary to insert boundary conditions which, in case of a static problem, represent specified forces and/or displacements at system edges.
Prior to transformation to the global coordinate system, the local stiffness matrix of an element can be expressed, in its explicit form, as: (3) (4) (5) Neno Torić, Alen Harapin, Ivica Boko It can clearly be seen from the above expression (8) that the local stiffness matrix, apart from depending on the beam length l , also depends on the material parameters: E, G, and the geometrical parameters: A, I y and I z .
In case of a real element (beam/column) subjected to a load, internal forces (and primarily bending and torsion moments) The material nonlinear analysis can very easily be applied by dividing the element into smaller parts (subelements), in which case the real stiffness is calculated for each subelement (Figure 1).

Numerical model for determining fire behaviour of structures
In order to conduct the nonlinear analysis, it is therefore necessary to define in which way the cross sectional stiffness at various levels of stress should be calculated.
2.3.Calculation of stress and strain and calculation of cross-sectional stiffness

Basic assumptions
Basic premises of the model for calculating stress-strain situation along the cross-section [12,18,19] are: -cross-sections remain flat even after deformation, -no sliding has been occurred at the connection of different materials, -uniaxial stress-deformation relationship (constitutive law) is known for all materials.

Cross-section deformation plane
The graphical view of a possible deformation plane, as compared to the previous state of equilibrium, is presented in Figure 2.An additional deformation of a given cross-sectional point is defined through equation of a plane: where: In the above expressions represents the vector of unknown parameters of the additional strain plane, and y, z the coordinates of the point in the Y-Z plane.The strain plane is described by its intersection with the X coordinate axis, and the components of relative rotations , around the Z and Y axis, respectively.If the considered section point has previous strain , its total strain is: The deformation is known and determined through equilibrium position via , by analogy to expression (9), by means of: If expressions (9) and ( 13) are inserted in the expression (12), we get: or: where is the vector of the total plane deformation parameters.

Relationship between stress and strain
The starting point is the known relationship between the uniaxial stress s and strain e for a given material, i.e. the so called constitutive law of material behaviour.For real materials, this relationship is basically curved, and is defined by uniaxial test or appropriate regulations.From the standpoint of numerical analysis, it is appropriate to define this relationship as linear by individual segments (Figure 3).The controlled error introduced in this way is negligible when compared to other assumptions.
The s-e relationship between any of the two i, j points of the diagram is described by the following expression: If the expression ( 14) is inserted in ( 16) and if the following change is made: then the stress in the sector under study can be described with the following expression: In the above expressions, the element E denotes the current modulus of elasticity of material (inclination of the straight line in the sector under study), while graphical interpretation can be seen in Figure 3.It should be noted that the stress is constant and defined for the known initial situation and for the assumed current deformation between the points i and j.The vector S v is composed of the vector of external forces S vp , which defines the initial deformation plane , and of the additional force vector , which causes the additional deformation plane .
To establish the balance, S v must be equal to S u , i.e. : In its developed form, the expression (24) is the system formed of three equations with an unknown vector 2.3.5.Bar reinforcement After determination of the total deformation of a given steel bar, the procedure continues by defining between which node deformations i , i+1 on the -diagram it is situated.Then the corresponding elastic model E, is defined, along with contribution of current mechanical properties of bar material, by analogy to expression (22), using: where A s denotes the area of the bar.All bars cross-sections are summed up.A part of the internal force vector the bar reinforcement is contributing to be defined by analogy with the expression (21), i.e. by: In case of bar-based materials, it is more appropriate to define the initial deformation with the discrete value , rather than with parameters of the initial deformation plane .

Larger area materials
The area of material with significantly greater area when compared to the area of the entire cross section, is assigned through convex polygonal elements without voids (finite Neno Torić, Alen Harapin, Ivica Boko In case of structures exposed to fire action, the stress and strain relationship of a material is highly dependent on temperature.For that reason, the behaviour of materials subjected to high temperatures is usually described as a set of stress -strain curves, each one for a specific temperature.At that, temperatures situated in between the specified ones are obtained by linear interpolation.

Equilibrium equation
The vector of cross-sectional internal resistance forces S u is a function of the resulting deformation plane and the s-e relationship of a specific material.If they are known, the S u can be calculated by integrating stress along the composite section area: where N u is the internal longitudinal force, M zu and M yu are the corresponding moment forces with respect to coordinate axes, W is the area of a specific material, and S summation across all materials m.If (18) is inserted in (19), then: where: Numerical model for determining fire behaviour of structures elements -FE).With the exception of bar reinforcement, only one type of material can be situated in the area occupied by a single FE.Each FE is defined with a list of nodal points and their coordinates, and with an index of material properties.In other words, contours of each material are first approximated with the polygon, and then the bounded area is divided into FEs (Figure 4).After determination of the resulting deformation plane on the observed FE and position of the corresponding neutral axis in the previous iteration, a set of straight lines parallel to it, with FE points along these lines, is defined, together with deformations equal to nodal deformations i of the working diagram.The points of intersection of these straight lines with sides of each FE is sought, and so the areas Ω ei (subelements) with constant modulus of elasticity E are defined on each of them.The matrix I e for each of these areas is shaped as follows: and it has been obtained by adding up all sides of subelement, using the expressions: where arranged pairs ( , ) and ( , ), represent coordinates of boundary points of the side under study, and n is the number of nodes (sides).The summing up is conducted across all sides of the subelement.A part of vector of internal forces, to which this subelement is attributed, according to expression (20), is defined with: where: Mechanical properties and a part of the internal force vector of one finite element are obtained by summing up appropriate properties of all areas on this element, and of individual materials through all finite elements that describe this material.Analogously, by summing across all materials we obtain total properties of the composite cross-section.The rise of temperature in space affected by real fire is dependent on a number of parameters: space covered by a combustible material, quantity of combustible material, size of ventilation openings, and dimensions of the affected space.Considering a great number of parameters influencing growth of temperature in space, when conducting experiments related to fire behaviour of structures, researchers the most often use the typified form of temperature growth in stove expressed/described by the fire growth curve defined by the following equation:
The spatial domain is approximated with a specified number of finite elements, as shown in Figure 5.The presented heat transfer model uses 8-node 3D finite elements.

Integration of discrete system equations
The system of non-linear ordinary differential equations ( 35) is solved through integration of equations between two neighbouring time points for a relatively small time interval .Temperatures at the beginning of the time interval are known and are used to calculate temperatures at the end of the time interval, i.e. at the moment t+ .By applying the mixed integration approach to the equation system (35), we obtain: where: -known temperature at the beginning of the time interval, -unknown temperature at the end of the time interval, -interpolation parameter, -field production vector at the end of the time interval, -boundary heat flow vector at the end of the time interval, -incremental time step.
The iterative procedure for determining the value starts by assuming the equivalence = , in the first iterative step within the time increment .After that, a new temperature value at the end of time interval is calculated, until the following condition is met: where: -incremental time step.

Transient heat transfer
The transient heat transfer is a time-dependant heat transfer process in which the temperature field created by transfer of heat within the area under study changes in time.The model implemented is based on the 3D transient nonlinear heat transfer model.The differential equation describing this process in a spatial setting is defined by the expression: where: -density of matter (kg/m 3 ) -specific heat capacity (J/kgK) -thermal conductivity coefficient tensor (W/mK) -heat diffusion tensor (m 2 /s) By applying the weak formulation to the equation ( 34), and by using the Galerkin method for selection of approximate solution function, we obtain the system of p ordinary differential equations: where: -domain, -edge of space domain, -shape functions for approximate solution, -total number of nodes in space domain discretisation.
The matrix C is called the capacitive matrix, matrix K the transmission matrix, vector F the thermal load vector, and Q the thermal inflow vector.The heat flow occurring as a result of fire is composed of convection and radiation.The heat flow on the surface of heated element is determined using the following expression: (42) Numerical model for determining fire behaviour of structures Figure 6 shows discretisation of a simple structure with the link between the spatial beam-element system, 2D and 3D network.The 2D network is used for discretisation of beam structure cross-sections, i.e. for calculating the stress-strain relationship along cross-section, and for stiffness calculation.The 3D network (model) is used for heat transfer calculation/analysis.The procedure starts by defining the spatial beam (frame) system, i.e. the initial and the final node must be defined for each beam (column/beam).In addition, the cross section with behaviour pattern ( -diagram) must also be defined for each beam element and for each material used (2D network).The 3D network for heat transfer calculation is then automatically generated along the beam element (beam/column).It is also necessary to define the number of subelements the final element will be divided into.
A cross section and the law of behaviour from the global element (beam/column) are attributed to each subelement.The heat flow calculation is made on the global 3D model and an instantaneous temperature is obtained in every node.The mean (average) temperature is calculated at every 2D network element (which represents the cross-section of the element), and the constitutive material law is corrected.After that the stress-strain situation in cross-section can be defined, as well as the cross-section stiffness, which represents stiffness of the beam (beam, column) subelement.
The model is incremental and linear along individual increments.
The procedure starts from the cross-section level (zero position) by calculating real cross-section stiffness values for unstressed cross-section, according to expressions ( 22), ( 25) and ( 27).The values obtained in this way represent the initial stiffness of crosssection.The initial stiffness values are used to calculate the initial beam-element stiffness matrix.It is significant to note that only axial and flexural properties of cross-section are set/ corrected by integration at the cross-section level, while shear and torsional properties are left unchanged.This is followed by the usual procedure of arranging the global stiffness matrix and the global load vector (2), and by system calculation.After calculation of internal forces at the end of the beam element, the position of the deformation plane is determined and new cross-section stiffness is defined.In general, two cases are possible: 1.The deformation plane position can be determined for cross section.In this case, the strength of the cross-section is sufficient to withstand the external force action (i.e. the action of forces obtained by linear analysis of the beam element system).2. The deformation plane position can not be determined for cross section, i.e. the procedure presented in section 2.2 is divergent.In this case, it can be stated that there is a failure at the cross-section, i.e. that local failure has occurred in the system, and hence that global system failure is also possible.
In concrete terms, the stiffness values are then made equal to zero, and the global analysis procedure is reiterated.
The procedure continues until the standard of the displacement increase vector does not fall under an arbitrary value, i.e.: In all practical cases, it can be assumed that the value amounts to 0.001.The incremental procedure and the computer program flow diagram are presented in Table 1 and Figure 6. (43) 1 The initial load vector, which is actually the nulvector (F=0) is assumed.On the basis of this load, the initial local stiffness matrices are determined for each beam element (20 and 22), and also the initial global stiffness matrix : This matrix represents the so called zero stiffness, i.e. the system stiffness not influenced by forces. 2 From the given external load on the beams/ columns, defining the forces on elements, and the vector of the applied forces: 3 Time loop, the first time step is set j=1. 4 The calculation of the temperature field in the 3D element (Fig. 5).For each cross section: calculation of the mean temperature and the correction of the material characteristics according to the calculated mean temperature.the convergence is satisfied, the procedure is finished and the results are printed.Procedure continues with new time step (4).If the convergence is not satisfied, the procedure continues-step (7).The value is an optionally chosen small value, usually 0,001.

8
The calculation of the new stiffness in 2D elements, according to corrected material characteristics and the corrected internal forces.
Procedure continues with new time step: j=j+1 and the calculation procedure continues to step (6).The procedure continues until certain accuracy is achieved or until divergence occurs.This divergence points out to that load capacity of the cross-section is reached, i.e. the failure of the element in that section and possible failure of the entire structure.An example of modelling experiment conducted by Wainman and Kirby [20] is considered in order to present capabilities of the described numerical model for predicting fire behaviour of structures.This experiment involves a freely supported steel element I 254/146, steel quality S275, length: 4.58 m, exposed to temperatures defined by the standard fire curve.A noncomposite concrete slab 12.5 cm in thickness and 65 cm in width, situated on top of the element, exerts load amounting to 2.21 kN/m along the length of the slab.Four concentrated forces K were taken as additional load.The forces amount to 32.5 kN and are placed in positions as defined in Figure 7, where we also see disposition of measurement points in which a time-related increase of temperature was monitored.The dimensions and geometrical properties of the steel section I 254/146 that were used as input data for modelling the freely supported steel element are presented in Table 2.
Numerical model for determining fire behaviour of structures    The comparison of results of vertical deflection measured in a half of the element, taken from the experiment conducted by Wainman and Kirby, with results obtained through the described model, is presented in Figure 14.

Figure 1 .
Figure 1 .Stress-strain and stiffness along the girder as related to the level of stress in cross-section

Figure 2 .
Figure 2. Graphical presentation of a possible deformation plane of the vector of internal forces which is obtained by integration of the stress along the entire area of the composite cross section.Matrix members I represent current mechanical properties of the cross section.

Figure
Figure 3. Possible stress -strain relationship for a material

Figure 4 .
Figure 4. Spatial discretisation of cross-section 2.4.1.Temperature rise in confined space and standard fire curve

Figure 5
Figure 5. a) Global discretization of the space frame.(b) The beam-column element.(c) The element's cross-section discretization.(d) The comparative body model for heat transfer analysis.(e) The stressstrain constitutive law of the element's cross-section.
of numerical model 3.1.Description of model and experiment

Figure 12 .
Figure 12.Comparison of temperature results obtained by experiment and model -top flange