Model parameter sensitivity and benchmarking of the explicit dynamic solver of LS-DYNA for structural analysis in case of ﬁ re

Due to the complex nature of structural response in ﬁ re, computational tools are often necessary for the safe design of structures under ﬁ re conditions. In recent years, use of the ﬁ nite element code LS-DYNA has grown considerably in research and industry for structural ﬁ re analysis, but there is no benchmarking of the code available in the ﬁ re science literature for such applications. Moreover, due to the quasi-static nature of structural response in ﬁ re, the majority of the computational structural ﬁ re studies in the literature are based on the use of static solvers. Thus, this paper aims at benchmarking the explicit dynamic solver of LS-DYNA for structural ﬁ re analysis against other static numerical codes and experiments. A parameter sensitivity study is carried out to study the e ﬀ ects of various numerical parameters on the convergence to quasi-static solutions. Four canonical problems that encompass a range of thermal and mechanical behaviours in ﬁ re are simulated. In addition, two di ﬀ erent modelling approaches of composite action between the concrete slab and the steel beams are investigated. In general, the results con ﬁ rm that when numerical parameters are carefully considered such as to not induce excessive inertia forces in the system, explicit dynamic analyses using LS-DYNA provide good predictions of the key variables of structural response during ﬁ re.


Introduction
Modern building designs and innovative architectural solutions pose a challenge to structural engineers. This is particularly the case for structural fire engineers due to the complex interactions of modern structural systems in fire. The performance of even generic structures exposed to fire is not straightforward and cannot be easily predicted. As a result, there is often the engineering need to be able to assess structural behaviour under fire conditions from first principles and not rely on blanket prescriptive guidance.
The behaviour of isolated structural elements under standard fire conditions through furnace testing has been extensively studied over the past decade, and can now be predicted with some degree of accuracy using analytical and computational means [1]. However, it has been shown in the past [1][2][3] that structural fire performance of isolated elements does not resemble the performance of a whole structure. The whole structure performance in a real fire depends on a number of factors. They include restraint, stress redistribution, composite action, and continuity within the structure [4]. The involvement of the many variables makes the analysis and prediction of the fire performance of realistic structures a difficult process. Standard fire tests provide unrealistic results [5]. They do not represent real fire conditions in the compartment and base fire resistance on the performance of the individual elements ignoring the effects of the surrounding structure. Conversely, full-scale testing of real structures is complex, expensive and time consuming. In addition, the limited number of full-scale experiments carried out worldwide (e.g. Cardington tests [5]) has been on buildings of generic rectangular geometry. Thus, they cannot be generalised to predict the performance of all structures, especially where more innovative irregular structural arrangements are used. As a results, designers use computational tools to predict and assess the performance of complex structures under fire conditions.
With increasing computational capabilities, the fire resistance assessment of various structural arrangements under different fire scenarios is becoming more and more used in practice. However, these models have to be benchmarked against experimental data or known solutions to make sure that they produce accurate and physically correct results. Most commonly used numerical models for structural fire analysis include commercial general finite element analysis packages (Ansys, Abaqus) and purpose-based finite element models developed or extended specifically for structural fire analysis (Vulcan, SAFIR and OpenSees). All of these models have been widely used for structural fire analysis in the recent years and have been validated against various small-scale and full-scale tests (e.g. Cardington) [6][7][8][9][10][11][12][13][14][15][16][17].
More recently, researchers and designers [18][19][20] have adopted LS-DYNA for structural fire analysis. Kilic and Selamet [19], and Selamet [18] used LS-DYNA to investigate the effect of fire location on the collapse of a 49 storey high-rise steel structure. The whole 3D building model was used for the analysis. Law et al. [20] studied the structural response of structural arrangements with bi-linear columns to fire. A slice of the 3D model of a generic substructure was used that included a 4-storey high bi-linear column with adjacent composite beams and concrete slabs. In all of these studies [18][19][20] the authors used beam and shell elements to represent steel beams and concrete slab respectively. However, neither of them included the benchmarking of the adopted LS-DYNA model. Both et al. [21] published a benchmark for predicting the structural fire response of a centrally loaded steelconcrete column. The column was modelled using solid elements. The results of the coupled thermal-mechanical analysis were presented for ANSYS, LS-DYNA, and ABAQUS. Temperature and displacement development in the column showed a good agreement between those software packages. The maximum differences between the peak values in relation to LS-DYNA results were approximately 23°C (5%) and 0.23 mm (23%), respectively. Kwaśniewski et al. [22] carried out a coupled thermal-mechanical analysis of a restrained steel column subjected to fire using LS-DYNA and validated against experimental results. The detailed 3D model (as in the benchmark by Both et al. [21]) adopted solid elements.
LS-DYNA is a commercial general purpose finite element software originally developed for highly nonlinear and transient dynamic analysis [23]. It is robust in the analysis of problems involving transient effects, contact and large deformations and has high computational efficiency. As a result, LS-DYNA is one of the most commonly used numerical explicit integration simulation programs. Common applications of LS-DYNA include automotive, aerospace, metalforming, and multi-physics problems. In structural engineering, common applications include earthquake, blast impact, and progressive collapse analysis. LS-DYNA has been used for the aircraft impact and progressive collapse analysis of the World Trade Centre (WTC) towers by NIST [24,25]. However, most likely due to the lack of available scientific work using LS-DYNA, for structural fire response analysis in the same WTC study a different numerical program was used. LS-DYNA as a software has been fully validated and verified by its developers (Livermore) for its generic applications. Even though, in recent years, the explicit dynamic solver of LS-DYNA has grown considerably in research and industry uses for the analysis of structures in fire [18][19][20][21][22], to the best of the authors knowledge, there is still no benchmarking of the code available in the literature specifically with regards to the structural fire performance. Available research is limited to the internal benchmarking work carried out in Arup [26] and detailed 3D solid based individual structural member models [22]. Solid elements are not frequently used for global models or for design purposes. In addition, due to the quasi-static nature of the structural response in fire, the majority of the structural fire analyses available in literature are carried out using static solvers.
Thus, this paper aims at benchmarking the explicit dynamic solver of LS-DYNA for the structural fire analysis of 3D composite structures and 2D steel frames against experiments and other numerical codes. An extensive parameter sensitivity study is carried out to study the effects of various modelling parameters on the kinetic energy and convergence to quasi-static solution. Four canonical problems that encompass a range of thermal and mechanical behaviours in fire are simulated. The term benchmarking is used in this paper to refer to verification and validation of computational models, based on the definitions adopted by the ASTM [27]. That is, evaluating the software for correct application to known benchmark problems and for physical correctness of the results [15,16]. The LS-DYNA model is benchmarked against fire tests on loaded steel framework results [28][29][30], two benchmarks published by Gillie [31] and results published by Rackauskaite and El-Rimawi [32] on the numerical study of 2D steel frames subject to localised fires.

LS-DYNA benchmarking models
For the benchmarking of LS-DYNA for structural fire analysis, we use the double precision LS-DYNA (Release 7.1.1) version. Each benchmarking case chosen for this paper encompasses different mechanisms of structural response in fire, which are required to get a realistic response. In total four benchmark cases are considered. The first benchmark is based on the natural fire test of the 2D steel frame carried out in 1987 [28][29][30]. It allows the benchmarking of LS-DYNA against experimental results and assessing whether the model correctly captures the effects of non-uniform heating, material non-linearity, restraint, and stress redistribution.
The remaining 3 benchmarks are based on and allow benchmarking of LS-DYNA against the results of numerical analysis published in the literature [31,32]. Gillie [31] has published results for 2 problems which provide benchmark solutions for structural fire analysis. They allow users to check whether their models capture the required phenomena, which occur when structures are heated. The first benchmark [31] is on a uniformly heated steel beam with 75% support stiffness. This benchmark allows to confirm whether the model captures the effects of material non-linearity, geometric non-linearity and restraint conditions [31]. The second benchmark [31] is on a heated composite concrete floor. This benchmark assesses whether phenomena such as stress redistribution, localised heating and composite action effects can be captured. These two benchmarks provide a computational challenge against most of the fundamental mechanisms that occur when structures are heated and, thus, were chosen for the benchmarking of LS-DYNA. In addition to the above, a third benchmark has been chosen based on the study by Rackauskaite and El-Rimawi [32] on the heating effects on a 2D steel frame. In the latter study non-uniform heating of the beams was assumed. Therefore, the adoption of this study as a benchmark allows to assess whether thermal bowing, restraint from the surrounding structure and stress redistribution are appropriately captured in the analysis. In the following sections these benchmarks are described and the results from the LS-DYNA analyses are presented.
2.1. Benchmark #1 (BM1): fire test on a loaded steel framework Benchmark #1 (BM1) represents a natural fire test on a 2D steel frame carried out in 1987 [28]. The test frame comprised of a 4.55 m long steel beam and two 2.53 m high columns. It was connected to secondary framework to prevent lateral instability. Details of the frame and loading are shown in Fig. 1. This test has already been successfully modelled in CEFICOSS [29], Abaqus and SAFIR [30]. Therefore, a similar modelling approach was used in LS-DYNA for comparative purposes.
Beam web temperature measurement from the test is only available at one instant during the whole fire exposure. Thus, member temperatures, which show a good agreement with the measured test data for beam flanges and column (see Fig. 1), were taken from [29]. The beam was subjected to non-uniform gas temperatures along its length. To account for this, a temperature reduction function following a sinusoidal shape was applied along the beam as in [29,30] with the temperature at the connection being 0.9 of the beam temperature at mid-span. Column temperatures along the height were assumed to be constant.
In LS-DYNA, the steel frame was modelled using Hughes-Liu beam elements with user defined cross-section integration. This element has a single integration point along its length in the middle of the element, where sectional plasticity is monitored. Hughes-Liu beam elements allow for the treatment of finite strains and are simple, computationally efficient and robust. It is the first beam element implemented in LS-DYNA and thus it is used as the default element type in the program [33]. The same beam formulation was used in [18][19][20]. A thermallysensitive steel material type MAT 202 formulation based on Eurocode 3 (EN 1993-1-2:2005) was used for the analysis. Thermal loading to the beam and columns was applied using the formulation (i.e. "load thermal variable beam") that allows the definition of variable through-thickness temperature distribution. Using this formulation, a known temperature-time history at different cross-section coordinates of beam elements can be defined. Along the length of the beam element only a single uniform temperature can be defined, due to the Hughes-Liu beam element formulation. Therefore, for this benchmark a uniform average temperature value over the beam element length was calculated and applied.
Due to symmetry, only half of the frame was modelled with the support conditions set as shown in Fig. 1. To represent the support from the secondary framework a discrete beam element was used at the column height of 3.2101 m with nonlinear elastic material (MAT 67) formulation. This material model allows defining 6 springs acting about the local degrees of freedom and nonlinear spring stiffness [23]. Spring behaviour was taken from [29,30]. LS-DYNA has a 3D model space and as a result, out of plane translational restraint to the nodes was added to simplify the benchmark into a two dimensional space.

Benchmark #2 (BM2): uniformly heated beam
Benchmark #2 (BM2) [31] represents a 1 m long restrained rectangular steel beam with cross-section dimensions of 35 mm x 35 mm. A uniform load of 4250 N/m is applied to the beam (see Fig. 2). The beam is linearly heated up to a temperature of 800°C and then cooled down to room temperature. The temperature distribution within the beam is assumed to be uniform. Different levels of restraint are considered and include 0% (i.e. simply supported), 5%, 10%, 25%, 50%, 75%, and 100% (i.e. pinned) of support stiffness.
In LS-DYNA, the steel beam was modelled using Hughes-Liu beam elements as for BM1. The elastic-plastic material (MAT 004) formulation was used to represent the steel. This material model allows the user to define temperature dependent material properties as provided in the original benchmark [31]. One end of the beam was set as pinned and on the other end a discrete beam element was used with a linear elastic material (MAT 66) formulation. This material model allows defining 6 springs acting about the local degrees of freedom [23]. The discrete beam element was used to represent varying support restraint conditions of the beam. To model the 75% support restraint, the translational stiffness in the discrete beam material model was set to 1.902×10 8 N/m (75% of the beam stiffness, EA/L) as suggested in [31]. Out of plane translational restraint was added to all nodes to ensure a two dimensional behaviour. Illustration of this benchmark (BM2) is shown in Fig. 2.

Benchmark #3 (BM3): composite steel-concrete floor
Between 1994 and 1996 a series of full-scale fire tests were conducted on an 8-storey steel frame within the UK Building Research Establishment at Cardington to gain a better understanding of global structural performance in fire [5]. Benchmark #3 (BM3) [31] represents a simplified idealisation of the first Cardington test. In this benchmark part of the 4.5 m x 4.5 m x 0.13 m composite steel-concrete floor is exposed to linearly increasing heating and cooling. The heated area and the geometry of the benchmark are shown in Fig. 3. The concrete floor slab is connected to primary and secondary steel beams. The heated part of the secondary beams has a uniform temperature distribution and reaches a peak temperature of 800°C. The heated part of the floor slab is assumed to have a linear through thickness thermal gradient with the absolute peak temperature values of 600°C at the bottom of the slab and 0°C at the top. The uniform gravity loading applied on the slab is 5.48 kN/m 2 .
For this benchmark two LS-DYNA models were developed to compare two different approaches used to represent the composite action between steel beams and the concrete slab. In Model A (BM3-A), the beam elements representing the steel beams and the shell elements representing the concrete slab were defined using separate nodes. To represent composite action, a tied contact between beam and shell , and comparison of beam temperatures at mid-span and column temperatures at mid-height measured during the test [28] with CEFICOSS numerical simulation results [29,30] that were used as an input in LS-DYNA analysis (right). elements was used. This contact uses a penalty based formulation and force and moment resultants are transferred by discrete spring elements [23]. In Model B (BM3-B), the beam and shell elements share the same nodes in order to account for composite action effects between steel beams and concrete slab. Beam and shell element formulations with offset option were used to offset the central axis of the elements to geometrically represent the position of the elements in relation to one another. The illustration of the geometry and the two models are shown in Fig. 3. Both modelling approaches assume that there is a perfect contact between the steel beams and the concrete slabs, and that shear stud failure does not occur during the fire. The remaining parameters and configuration (e.g. mesh size, material properties, loading) were the same in both Model A and Model B. For primary and secondary steel beams, the Hughes-Liu beam element formulation was used. The concrete slab was modelled using the Belytschko-Tsay shell element formulation, which is computationally efficient. It is based on a combined co-rotational and velocitystrain formulation [33]. Material type MAT 172 and MAT 202 formulations were used for concrete and steel respectively. The steel material model is based on Eurocode 3 (EN 1993-1-2:2005) and the concrete model is based on Eurocode 2 (EN 1992-1-2:2004). Both material models are thermally-sensitive and allow the users to define temperature dependent properties by overriding the values from the Eurocodes. MAT 172 can be used to represent plain concrete, plain steel or a smeared combination of concrete and steel [23]. In this benchmarking study two materials were created using this formulation for plain concrete and plain reinforcement. These were then used to represent different through-thickness layers of the concrete slab section. The slab was divided into a chosen number of layers, each representing either plain concrete or plain reinforcement, with the middle layer characterising steel reinforcement. The equivalent thickness of the reinforcement layer (0.283 mm) was calculated based on the total volume of steel rebars in the concrete slab.
The thermal load to the beams and shell elements was applied using the formulations that allow the definition of varying through-thickness temperatures. It was identified by Gillie [31] that gravity loading was included in the uniformly distributed load on the concrete slab and thus it was not explicitly applied in the current model. The concrete slab was modelled as continuous with appropriate lateral and rotational constraints around the boundary.

Benchmark #4 (BM4): 2D steel frame
The fourth LS-DYNA benchmark model used for benchmarking, Benchmark #4 (BM4), is based on the numerical study by Rackauskaite and El-Rimawi [32] on the structural response of 2D frames subject to localised fires. The latter study [32] was conducted using an in-house specialist finite element software based on the secant stiffness approach that was developed specifically for structural fire modelling purposes [32,[34][35][36]. A three-storey three-bay steel frame was subject to six  Fig. 4. The heated columns were assumed to have uniform temperature distribution. The steel beams were subjected to a temperature gradient through the section. For the top flange, web and bottom flange of the heated beams temperature scaling ratios of 0.85, 1, and 0.86 were used respectively. This is to account for the increased temperatures of the web compared to the bottom flange due to its thickness and the heat sink effect of the slab. The structural members were linearly heated until the failure of the frame occurred.
As for the previous benchmarks, the steel beams were modelled using the Hughes-Liu beam element formulation. The element lengths for beams and columns were 1 m and 0.875 m respectively. This mesh density was used in [32] and thus the same mesh density was adopted in the current study as well for comparative purposes. Steel material type MAT 202 formulation was used for steel beams and columns with the default Eurocode temperature dependent material properties. Thermal loading to the beams was applied using the formulation that allows the definition of variable through-thickness temperature distribution. The uniformly distributed load on beams of 41.1 kN/m was converted to equivalent nodal loads in accordance to the original study [32]. Gravity loads (i.e. self-weight of the beam) are included in the latter load and, therefore, were not considered separately. Column base boundary condition was set as fixed. Out of plane translational restraint was added to all nodes to ensure a two dimensional behaviour.

Parameter sensitivity
Prior to benchmarking of LS-DYNA for structural fire analysis, a model parameter sensitivity study was carried out in order to determine the optimum balance between accuracy and computational efficiency. The effect of parameters such as load time scaling, mesh density, number of beam and shell element integration points on the convergence of the solution was investigated. The initial model parameters are identified in Table 1 and are described in the following sections. During the parameter sensitivity study only one parameter was varied at a time. The other model parameters were kept constant as shown in Table 1. Based on this parameter sensitivity study, the final model parameters were then chosen for benchmarking and the comparison of the LS-DYNA analyses results with the test or model outputs as provided in [30][31][32].

Time scaling
The response of structures to fire is frequently idealised as a quasistatic problem [37]. This is due to the relatively slow increase in temperatures and therefore the rate of change of displacements (until a failure occurs) of structural components during a fire. Dynamic instabilities generated during fire are relatively small compared to the natural frequency of a structure [37]. In addition, commonly used material models for structures in fire are temperature dependent only with creep considered implicitly. Therefore, the heating or cooling rate does not influence the modelled material response. LS-DYNA is a nonlinear transient dynamic finite element code. Thus, it uses real-time to solve the equation of motion (Eq. (1)) [33].
where m is the mass, ü is the acceleration, c is the damping coefficient, ui s the velocity, f u ( ) int are the internal forces, and p t ( ) are the external forces.
The modelling of structural response in fire using the real time fire duration would require a significant computational time. Due to the mainly quasi-static nature of the structural response to fire problem, when a dynamic analysis is carried out, the most common approach is to scale-down the "real" time in the model and therefore apply static and thermal loads over a shorter period. However, scaling the load application time implies that the scaled down "load" will be applied faster. Thus, for the LS-DYNA explicit dynamic solver (or any other dynamic solver), it could lead to the introduction of significant inertia forces to the system, which could impact the structural response. During the thermal load application inertia effects could be initiated by thermal expansion and large deflections due to high temperatures developing over a short period of time.
A general approach to avoid inertia effects and, thus, high dynamic oscillations is to ramp the load linearly over a period of time to minimise any numerical instabilities and to reach a near to steady-state solution (due to static loads) or quasi-static solution (due to thermal loads). This approach can be combined with the application of damping, which dissipates the kinetic energy in the system and allows the use of shorter preload and thermal load durations and thus computational time. Other approaches could also involve mass scaling, however, this approach is not assessed in this work. In this paper we investigate the effect of preload and thermal load application duration (see Fig. 5) and damping on the kinetic energy and force development in LS-DYNA for the four benchmark cases identified in the previous section. Kinetic energy gives an indication on the amount of inertia force present in the system. A general rule of thumb is that a solution is considered to be steady state and quasi-static when the kinetic energy to internal energy ratio in the system for most of the time is less than 0.1% and 5% [38][39][40][41], respectively. In this study the same criterion is used based on the 95th percentile of kinetic to internal energy ratio during the fire exposure.
Preload durations, t pre , are varied between 1 and 64 s for BM1, BM2, and BM4, and between 1 and 256 s for BM3. The thermal load application duration, t h , up to a maximum temperature is varied from 0.5 to 512 s for all benchmarks. The upper boundaries used are based  on the convergence of the output and kinetic to internal energy ratios for different models. The preload is applied first and then kept constant for 1 s. Following that once the steady-state solution has been reached, the thermal load is applied. In the cases where cooling is considering, cooling duration to ambient temperature, t c , is the same as the heating duration, i.e. t t = c h . During the thermal load application period, the static load is unchanged and considered to be constant. In the cases where damping is applied during preload only, the critical damping factor is reduced to 0 before the application of the thermal load. The order of load application is shown schematically in Fig. 5.
In addition to varying the preload and thermal load duration, the application of global damping was considered. For each benchmark, an LS-DYNA implicit eigenvalue analysis was run to determine the lowest fundamental frequency of the structural system. Material type MAT 202 formulation is not currently available on the implicit solver of LS-DYNA R7.1. Thus, to carry out the eigenvalue analysis for BM1, BM3 and BM4 an elastic material (MAT 1) formulation for steel parts was used. The resulting frequencies, f , were 62.4, 77.06, 7.70, 7.71, and 4.28 Hz for BM1, BM2, BM3-A, BM3-B, and BM4, respectively. These frequencies were used to determine the critical damping factor, d cr , for the structural system using Eq. (2) as suggested in [23]: The static load was applied using a half sine function shape over the duration 2 times the normal period (1/f ) of the system in the damped cases. After preload, the load was kept constant for 1 s. It should be noted that there may be cases where a dynamic analysis would be best suited for structural fire analyses such as in the case of fire induced progressive collapse [37,42] and therefore selecting the appropriate analysis parameters that would minimise the introduction of pseudoinertia effects is crucial.

Mesh density
One of the key model parameters that can influence the accuracy of the numerical solution is mesh density. Thus, for all the benchmark cases, a mesh convergence study is carried out. For cases BM1, BM2, BM3, and BM4, the beam element length is varied between 0.03 and 0.476 m, 0.00625 and 0.5 m, 0.0625 and 0.5 m, and 0.0625 and 4 m, respectively. Column element length for cases BM1 and BM4 is varied between 0.048 and 0.386 m and between 0.0625 and 0.875 m, respectively. The beam and shell element lengths for case BM3 are the same. For the mesh density study of BM1, a uniform temperature along the beam length is assumed. This is to avoid any secondary effects due to the local increase of temperatures with a finer mesh size. Hughes-Liu beam elements have one integration point in the middle of the element and therefore a uniform average temperature value has been assumed along the length of each element as identified previously in Section 2.1. With finer mesh, i.e. a higher number of elements with a shorter beam element length, the average temperature over the element length applied to beam elements would increase in case of non-uniform heating along the length of the member. This would therefore have an effect on the thermal expansion and stresses within the member. Thus, the outcome of the sensitivity study would not be only mesh dependent but also temperature dependent. Uniform heating was only considered as part of the mesh sensitivity study.

Beam element integration
As noted earlier, in LS-DYNA, the Hughes-Liu beam element formulation has a single integration point in the middle of the element and user defined cross-section integration [33]. Integration points for a rectangular cross-section (BM2) and an I-section (BM1, BM3, and BM4) are defined as shown in Fig. 6. The number of integration points can vary depending on the desired accuracy required on the throughdepth plasticity. Each fibre has a prescribed uniaxial stress-strain relationship. A greater number of integration points can also more accurately represent the thermal gradients within a structural member. In this study, the beam cross-section integration refinement factor, k (see Fig. 6), for all cases was varied between 0 and 11. That is, the number of cross-sectional integration points was varied for the beam elements between 9 and 169 for the rectangular cross-section (BM2), and between 6 and 50 for the I-section (BM1, BM2, and BM3).

Shell element integration
In LS-DYNA, based on the modelling approach presented for BM3 in this paper (see Section 2.3), the concrete slab is modelled as a composite section. Different shell through-thickness integration points are defined for steel rebars and slab with appropriate material models. A smeared cracking approach is adopted and therefore rebar rupture is not explicitly considered and would need to be assessed by the modeller independently. The rebar layer is modelled at the centre of the shell Fig. 5. The sequence of load application used in the LS-DYNA modelling approach for structural response to fire analysis. t pre is the preload duration, t ss is the steady-state, t h is the heating duration, and t c is the cooling duration.  Fig. 7). The total number of layers defines the number of integration points. In this study, the number of shell integration layers was varied between 3 and 11. The upper boundary used was based on the convergence of the output. In all cases only one layer at half thickness of the shell element was used to represent the steel rebars. In the original benchmark, [31] the exact location of the rebar in the concrete slab is not defined. The rebar position used in this paper was chosen based on the previous work by Gillie [43] on Cardington tests that this benchmark originated from and a study where the geometry of BM2 has been used as one of the investigated cases [44]. In the latter work, the depth of rebar in a concrete slab with trapezoidal decking was indicated as 55 mm, which is approximately at the centre of the concrete slab.

Effect of preload due to gravity and static load
The variation of the kinetic to internal energy ratio, displacement and axial force with preload duration during the steady-state (t ss ) are shown in Fig. 8 for all benchmark cases. It can be seen that with an increased preload duration the amount of kinetic energy reduces. The same load is applied over a longer period and thus the pseudo-inertia effects which are associated with the time scaling decrease. For BM1, BM2 and BM4, the kinetic energy present in the system is very low, in the range of 1×10 −4 % and 2×10 −3 % of total energy in the system even for the shortest preload duration of 1 s. That is less than 0.1% and indicates a steady-state solution. Some oscillations can be seen in the heated beam mid-span displacement and axial force results for preload application durations up to 32 s in BM1, BM2 and BM4. However, their amplitudes are very lowup to 0.01 mm (BM1), 0.007 mm (BM2) and 0.06 mm (BM4) for displacements, and up to 115 N (BM1), 10 N (BM2) and 150 N (BM4) for axial forces for a preload time of 1 s. Thus, the impact of the preload duration (t pre ) can be considered to be negligible for these benchmarks.
On the other hand, for BM3-A and BM3-B, the results show a much higher sensitivity of the heated beam mid-span displacements and axial forces with varying preload application duration. The kinetic energy is below 5% of the internal energy for the two cases and is decreasing as expected with increasing preload duration. A steady-state solution (kinetic / internal energy < 0.1%) is only reached for preload times of 64 s (BM3-A) and 128 s (BM3-B). The resultant beam mid-span displacements and axial forces show oscillations and average values changing with time even after preload. Unlike in BM1, BM2 and BM4, for BM3 the displacements and axial forces do not oscillate around a similar steady-state solution for varying preload durations. This may be due to the less confined boundary conditions in the model as vertical restraint is only provided at the centre of the concrete slab with continuity restraint conditions around the edges and due to the larger mass. Large slab mass and displacement rate contribute to higher inertia effects. In addition, the steady-state displacement solution for the highest preload time of 64 s is similar for both cases BM3-A and BM3-B and is around 3.1 mm. However, the difference in the axial forces is quite high, approx. 10 kN, considering that both cases represent the same benchmark problem.
When damping is applied, the kinetic energy in all cases is close to 0%. Results do not indicate any oscillations and show a nearly perfect steady-state solution. This solution in cases BM1, BM2 and BM4 is almost identical to results where damping was not applied. Nevertheless, this is not the case for BM3. The damping solution differs by up to 5 kN and 2 kN when compared to the solution with no damping for cases BM3-A and BM3-B, respectively.

Effect of thermal load
The variation of the relative error, 95th percentile of the kinetic energy ratio, displacement and axial force with preload and thermal load application duration are shown in Fig. 9. The relative error is based on the results of the case with longest preload (i.e. 32 s) and longest thermal load (i.e. 512 s) application duration for each benchmark. This is because for this case the inertia effects are the smallest and the solution is expected to be the one that is closest to quasi-static. It can be seen from Fig. 9 that thermal load duration induces larger inertia forces in the system compared to preload duration (see Fig. 8).
For the same application duration of static and thermal load of 1 s kinetic energy excited in the system for the latter is between 0.3% and 450% for all cases, while for static load it is only between 1×10 −5 % and 5% of the internal energy. The quasi-static solution based on the kinetic to internal energy ratio criteria ( < 5%) is reached for the time scale factor of 100 (BM1) and heating durations of 0.5 s (BM2), and 1 s (BM4) (see Table 2). However, axial forces for cases BM1, BM2 and BM4 indicate dynamic oscillations even though kinetic to internal energy ratio is below 5% for a time scale factor of 1000 (BM1) and a thermal load duration of 1 s (BM2 and BM4). Based on the relative error and convergence of the results a quasi-static solution is only reached for a time scale factor of 100 (BM1) and heating duration of 16 s (BM2 and BM4). At these values kinetic energy ratio is less than 0.01%. This indicates that the generally applied ad-hoc rule that the kinetic to internal energy ratio needs to generally be less than 5% to achieve a quasi-static solution might be too high for structural analysis in case of fire although further research would be required to make generalised conclusions.
For BM3-A and BM3-B, the peak mid-span displacement and axial force do not show a clear convergence. Also, some influence of the preload duration can be seen, unlike in BM2 and BM4. Yet, this influence seems to be smaller in comparison to the effect of the thermal duration. An interesting observation is that for case BM3-A the kinetic to internal energy ratio does not decrease with higher heating duration (between 2 and 32 s) as in other cases. For BM3-A, the peak beam midspan displacement is converging to approximately 0.06 m higher value in comparison to the results of BM3-B. This indicates a stiffer response for BM3-B due to using shared nodes for beam and shell elements. For the cases of BM3-A and BM3-B the quasi-static solution based on the kinetic energy criteria is only reached for a thermal load application duration of 128 s.

Mesh density convergence
The results for all cases assessed are shown in Fig. 10. Relative errors were calculated based on the difference to the finest mesh. It can be seen from Fig. 10 that for cases BM1, BM2 and BM4, the results converge with decreasing element size as it would be expected for a correct solution. Result convergence is reached for beam element sizes of 0.24 m, 0.05 m and 0.5 m for cases BM1, BM2 and BM4, respectively, with relative error less than 2%. However, for case BM3 as for the time scaling study, results do not show clear convergence. In addition, for the finest mesh size, oscillations and instabilities (BM3−0.0625 m) and an early failure (BM1−0.06 m) in the axial force development can be observed. This is likely due to the development of high localised stresses because Hughes-Liu beam elements have only one integration point along the length. Therefore, with a decreasing element size the constant stress along the length would increase. There is no clear relationship between the kinetic energy in the system and the change in element length comparing the different cases. Due to the instabilities present in the solutions for the finest mesh densities, the mesh with the beam element length of 0.25 m is selected as the "converged" solution for BM3. It should be noted that these analyses are based on capturing the global structural fire behaviour. The selected mesh densities for these problems would likely need to be different if localised responses such as concrete cracking are desired to be captured.

Effect of the number of beam element integration points
The variation of the relative errors of the beam peak mid-span displacement and axial force, kinetic energy and heated beam mid-span displacement development with the change of the integration refinement factor is shown in Fig. 11. The errors were calculated based on the numerical results for an integration refinement factor of 11, which includes the maximum number of fibres. For all cases except for case BM3-A, a convergence of the resultant values with an increasing integration factor can be observed. Yet, the relative errors for BM1, BM3-A, BM3-B and BM4 are negligible ( < 3.5%) as it can be seen from the development of beam mid-span displacement (Fig. 11). This indicates that increasing the number of beam integration points for these models (BM1, BM3, and BM4) does not improve the solution significantly. However, this is not the case for BM2. The development of beam mid-span displacement and the resultant values are affected by beam integration refinement. Convergence is only reached for a refinement factor, k, of 5 which gives a more refined solution in comparison to the case when k is 0. In the cases with k=1 and k=3 relative error reaches up to 20%. Thus, in some cases the integration factor has to be carefully chosen in order to obtain conservative results.
The kinetic energy in all cases shows negligible variation with the change in the number of beam integration points. The maximum difference is less than 1.5%. This indicates that the number of beam cross-section integration points is unlikely to affect the inertia of the system. It should be noted that if there is a scenario such that the response is significantly underestimated with a chosen number of integration points, this assumption may need to be reassessed.

Effect of the number of shell element integration layers
The variation of the key structural variables with the number of shell element integration layers is shown in Fig. 12. The relative errors, as in the previous sections, have been calculated based on the outcome for the model with the largest number of shell element through- thickness integration layers (i.e. 11). The results aashow a convergence with an increasing number of slab layers. The difference in the peak axial forces between the models with lowest and highest number of slab layers is up to 75% and 39% for BM3-A and BM3-B, respectively. However, these high values of error may be misleading. The cases with the number of shell element through-thickness integration layers between 5 and 11 have significantly higher values of kinetic energy present in the system compared to the models with only 3 integration layers. The differences in kinetic energy are between 12.5 to 16% (BM3-A) and 6-12% (BM3-B) of the internal energy. This indicates sufficiently larger inertia that could influence the results. On the other Fig. 9. Variation of relative error of peak displacement and axial force, 95th percentile of kinetic to internal energy ratio, displacement (preload 32 s), and axial force (preload 32 s) with preload and thermal load application duration for all benchmark cases. Plotted data is of the following structural members: BM2mid-span displacement and axial force; BM3heated beam mid-span displacement and axial force; and BM4mid-span displacement and axial force of the ground floor edge beam. DH refers to the cases where damping was continued during the heating. Table 2 Quasi-static solution details based on the kinetic to internal energy ratio criteria ( < 5%). hand, the axial force development for the 3 integration layers indicates some instabilities between temperatures of 300 and 400°C. Thus, for the purposes of this study the convergence is assumed to be reached for 7 shell element integration layers.

Benchmarking of LS-DYNA
In this section, the final results for each benchmark case obtained with LS-DYNA are presented and compared with the results published by Cooke and Latham [28], Santiago et al. [30], Gillie [31], and Rackauskaite and El-Rimawi [32]. Based on the time scaling and parameter sensitivity studies, the final model parameters were chosen and are shown in Table 3. The parameters were determined based on the convergence of results and computational time efficiency. The results for BM1, BM2, BM3, and BM4 are shown in Fig. 13, Fig. 14,  Fig. 17, and Fig. 18, respectively. The LS-DYNA final solution errors in respect to other software packages and/or experimental results are summarised in Table 4.
Benchmark BM1 represents the experiment on the 2D steel frame heated until failure. The results shown in Fig. 13 show a good agreement between the LS-DYNA solution and the experimental results [28] for the beam mid-span deflection and deflected shape of the frame at 16 min. Predictions from LS-DYNA on the deflections also compare well with other software packages (CEFICOSS [29], SAFIR and Abaqus [30]) and indicate a conservative solution. However, in comparison to the latter software packages, LS-DYNA underestimates peak axial forces by approximately 30 kN (25%). The sensitivity study presented in the previous sections indicates that there is very little or no influence from the investigated parameters on the development of axial forces. Therefore, LS-DYNA errors in comparison with other software packages could be associated with the inherent assumptions in the programs and the element types used to represent the support from the secondary framework. In the case without the lateral support (i.e. that discrete elements are not used), LS-DYNA underestimates peak axial forces by only approximately 2.9-4.3 kN (6.7-11.7%) in comparison with CEFICOSS, SAFIR and Abaqus [30]. It should be noted, however, that LS-DYNA gives a conservative prediction of the failure time by up to 1.4 min.
In benchmark BM2 a single rectangular steel beam with 75% restraint was heated and then cooled. For the BM2 model as seen in Fig. 14 LS-DYNA predicts well the development of heated beam midspan displacements and axial forces with temperature in comparison to other software packages (Abaqus, Vulcan, and Ansys [31]). Unlike in BM1, in BM2 the peak axial force predicted with LS-DYNA is within Fig. 10. Variation of relative error of peak displacement and axial force; 95th percentile of kinetic to internal energy ratio; and axial force with mesh size for all benchmark cases. Plotted data is of the following structural members: BM2mid-span displacement and axial force; BM3heated beam mid-span displacement and axial force; and BM4mid-span displacement and axial force of the ground floor edge beam.  Fig. 15. The results indicate that with the final chosen model parameters, LS-DYNA is able to predict well the development of deflections and axial forces for all restraint levels except for the 0-10%. For the latter cases high dynamic oscillations occur indicating large inertia effects for very low restraint levels. An interesting observation is that for these cases inertia effects could only be dissipated to achieve a quasi-static solution as in Abaqus by applying damping during the heating duration (see Fig. 15). The comparison of the 95th percentile of the kinetic to internal energy ratios (see Fig. 16) shows a significant reduction (a factor up to 30,000) in kinetic energy in the system when damping is applied for cases with low restraint levels (0-10%). Results shown in Fig. 16 also indicate that the application of damping during the heating duration may not always be the most effective or efficient solution. For restraint levels higher than 10% solutions with damping have very similar level of kinetic energy in the system as solutions without damping. The maximum difference in kinetic to internal energy ratio is up to 0.00006%. In benchmark BM3, part of the composite concrete floor slab and the supporting beams was heated and then cooled. For BM3, the LS-DYNA prediction in comparison to the prediction by Abaqus [31] is worse than for BM1 and BM2 (see Fig. 17). While LS-DYNA is able to predict the development of axial forces during the heating for BM3, it underestimates the values during cooling below the temperature of 400°C. LS-DYNA is unable to capture the peak axial force of approximately 1.5 MN at the end of cooling. The predicted values are 293 kN (19%) and 394 kN (26%) lower for BM3-A and BM3-B, respectively. Interestingly, this value has been captured by the model with the lower value of shell through-thickness integration points (i.e. 3) (see Fig. 12). Comparing models BM3-A and BM3-B, the latter one gives a closer prediction of axial force development in the beam during heating in comparison to Abaqus. In BM3-B a shared node between beam and shell elements was used to represent the composite action. Both BM3 models (A and B) show significantly higher beam mid-span deflection development values than Abaqus. Displacement results using BM3-A and BM3-B models are by up to approximately 1.4 and 1.2 times higher, respectively, compared to Gillie's [31] results. Thus, in the terms of deflections LS-DYNA solution can be considered as conservative.
The differences in LS-DYNA and Abaqus solutions could be due to different inherent assumptions in the analysis programs or input parameters used. No details of the input file for Abaqus [31] for BM3 Fig. 11. Variation of relative error of peak displacement and axial force (N), 95th percentile of kinetic to internal energy ratio, and displacement development with beam integration refinement factor for all benchmark cases. Plotted data is of the following structural members: BM2mid-span displacement; BM3heated beam mid-span displacement; and BM4mid-span displacement of the ground floor edge beam.   [32]. Fig. 13. Comparison of the experimental BM1 results [28] and finite element analysis results from CEFICOSS, SAFIR and Abaqus [29,30] with results obtained using LS-DYNA. Variation of heated beam mid-span deflection (left) and axial forces (middle) with temperature and the deflected shape of the frame at 16 min. Deflected shape scale factor is 10. have been published. Moreover, some details on the benchmark were not clearly defined in [31], for example, the location of steel rebar in the concrete slab and whether the self-weight of the members is to be included. Considering the uncertainty in the benchmark itself, Abaqus parameters used for the analysis, and significant sensitivity of the LS-DYNA solution to various parameters investigated in this paper, it is difficult to identify more specific reasons for the differences in LS-DYNA and Abaqus solutions. In addition, LS-DYNA models BM3-A and BM3-B at the end of heating give considerably different deflected shapes, especially in the bays adjacent to the heated slab. Deflections in these unheated bays in BM3-B are almost as high as in the heated part of the slab and thus they do not seem to be realistic. Gillie [31] has not published the deflected shape of the whole model for benchmarking. In benchmark BM4 various bays of the 3 storey × 3 bay 2D steel frame [32] were subjected to heating. The comparison of the results from [32] and the LS-DYNA solution for the different cases is shown in Fig. 18. For most localised fire cases the results show a good agreement. The development of the beam mid-span displacements and bending moments is almost identical in LS-DYNA and [32]. For axial forces LS-DYNA produces slightly higher axial forces (up to approximately 15 kN or 15% for cases A to D), but they can be considered to be conservative. Also, the trend of axial force development is the same in the results produced by LS-DYNA and [32]. Nevertheless, for the BM4 cases where the top floor and middle floors were heated, even though LS-DYNA  predicts similar general trends as in [32], the values predicted after a temperature of approximately 350°C are considerably higher than in [32]. At the temperature of 500°C the difference in the displacement, axial force, and bending moment values is up to 2.94 mm (24%), 47.3 kN (16%), and 16.1 kNm (9%), respectively. It is difficult to judge which model gives a better prediction but the results from LS-DYNA seem to be more conservative. Considering all four cases, LS-DYNA is shown to be capable of predicting the general behaviour of structures subjected to fire and to capture key phenomena such as the effects of material and geometric non-linearity, restraint conditions and stress redistribution. The maximum differences between LS-DYNA and other numerical analysis programs are up to 28% (42% for BM3-A) and 25%, for the prediction of deflections and axial forces, respectively (see Table 4). The highest differences are observed for the cases where symmetry or continuity boundary conditions were used (i.e. BM1 and BM3). In general, any discrepancies are likely to be of numeric nature due to the differences in finite element packages and their algorithms and are within acceptable limits. During cooling predictions are slightly worse and thus, LS-DYNA should be used cautiously when cooling is considered until further research is carried out. In some of the cases in order to achieve stable results and eliminate the effects of inertia long scaled heating durations had to be used. This led to very long computational times, which might be impractical sometimes. This could be overcome by using mass scaling. However, mass scaling should be used carefully. If it is too large additional inertia might be excited, which might have a significant influence on the final result.

Conclusions
In recent years, LS-DYNA has grown considerably in research and industry uses for the analysis of whole structures in fire. However, a benchmarking of the explicit dynamic solver of LS-DYNA for structural fire analysis has not been published in literature. Thus, in this paper, a benchmarking study of the explicit dynamic solver of LS-DYNA for the analysis of the structural response in fire and a model parameter sensitivity and convergence to quasi-static solution studies have been carried out. Four different canonical problems have been modelled in LS-DYNA based on various experimental and numerical studies published in the literature.
The results of this study indicate that the explicit dynamic solver of LS-DYNA is able to capture the key phenomena of heated structures. For all the benchmarks, it is able to predict the development trends of displacements, axial forces, and bending moments with increasing temperature within acceptable level of accuracy. However, during cooling in the composite concrete slab model, the LS-DYNA solution differs from the published results below the temperature of 400°C and is less conservative in predicting tensile axial forces. It is unknown whether this is because of the inherent assumptions in LS-DYNA or whether it is related to different input used. The predictions for the rectangular steel beam (BM1) during cooling compare well with other software packages which implies that LS-DYNA correctly captures the unloading behaviour of the steel material during cooling.
The time scaling sensitivity study indicated that inertia effects generated during the application of the thermal load to the system have  [31] with results obtained using LS-DYNA for varying levels of support restraint stiffness. Variation of heated beam mid-span displacement and axial forces with temperature for the models without (left) and with (right) damping included during heating. a more significant influence on the final result of the LS-DYNA explicit solution than preloading due to static loads. The high variation in the results and kinetic energy with different parameters, especially for the composite concrete floor slab case, highlights that an extensive parameter sensitivity study has to be carried out for every model to ensure that the LS-DYNA solution converges and is quasi-static. In addition, the time scaling study shows that the generally applied ad-hoc rule in explicit dynamic analysis that the kinetic to internal energy ratio for the most of the time has to be less than 5% to achieve a quasi-static solution might be too high for structural analysis in case of fire.
In this study LS-DYNA (Release 7.1.1) version was used. It should be noted that results obtained using a different LS-DYNA version might Fig. 17. Comparison of BM3 results published by Gillie [31] with results obtained using LS-DYNA models BM3-A and BM3-B. Variation of the heated beam mid-span displacement and axial forces with temperature (top), deflected shape at the end of heating from LS-DYNA model (bottom). Displacement scale factor is 5. Fig. 18. Comparison of BM4 results published by Rackauskaite and El-Rimawi [32] with results obtained using LS-DYNA. Variation of heated beam mid-span displacement (left), axial forces (middle) and bending moments (right) with temperature. For the references of beams in the frame refer to [32]. be to some extent different. Thus, it should be used cautiously and the model should be benchmarked before carrying out the analysis of interest. Additionally, further research is required to assess other element types such as solid elements or formulations of beams or shell elements different to those examined in this paper.