Recommended Finite Element Formulations for the Analysis of Off- shore Blast Walls in an Explosion

This study suggests relevant finite element (FE) formulations for the structural analysis of offshore blast walls subjected to blast loadings due to hydrocarbon explosions. The present blast wall model adopted from HSE (2003) consists of a corrugated panel and supporting members, and was modelled with shell, thick-shell, and solid element combinations in LS-DYNA, an explicit finite element analysis (FEA) solver. Stainless and mild steels were employed as materials for the blast wall model, with consideration of strain rate effect throughout ten (10) pulse pressure load regimes. The obtained FEA results were validated by experimental data from HSE (2003) with decent agreement. In the present study, recommended FE formulations with additional hourglass control functions were widely discussed from the perspectives of solution accuracy and computational cost based on a statistical approach. The obtained outcomes could be used for the structural analysis and design of offshore blast walls in the estimations of maximum and permanent deformations under blast loadings.


INTRODUCTION
Hydrocarbon explosion in offshore oil and gas installations is a highly unfavourable event that raises public concerns regarding safe practice in the industry. Among the most significant historical incidents are Piper Alpha in 1988, which sacrificed the lives of 167 crew members and Deepwater Horizon in 2010, during which approximately 4.9 million barrels of oil spilled into the Gulf of Mexico (GOM) resulting in a long-term environmental damage. Most recently, new procedures for structural assessment of offshore structures damaged by fire and explosion have been proposed by Kim (2014) and Kim (2016), in respectively.
Of several existing practical measures, corrugated blast walls are commonly installed as an integral part of offshore topsides as passive protection barriers to isolate hazardous hydrocarbon handling modules from personnel and critical equipment on board. As such, the structural integrity of these walls is critically emphasized throughout the design life of the entire platform structure. Hitherto, there are no unanimous guidelines for the design of explosion-resistant blast walls, though Technical Note 5 issued by the Fire and Blast Information Group (FABIG) has generally been referred to in designing stainless steel corrugated blast walls (FABIG, 1999). Within the industry, numerical methods like finite element analysis (FEA) are widely adopted over experimental and analytical methods, given its capability in modelling problems involving high structural complexities with various loadings and boundary conditions, in addition to providing greater insights into failure progression, making it the most cost-effective design and analysis tool.
With regards to blast wall structural assessment, significant research findings were presented by Langdon and Schleyer (2005a;2005b;2006) and HSE (2000) on blast response assessments of reduced and full scale stainless steel corrugated blast walls, through experimentation, analytical, and numerical studies. Louca et al. (2004) summarized the advantages and limitations of analytical single degree of freedom (SDOF) or the Biggs' method (Biggs, 1964) and numerical nonlinear finite element method (NLFEM) for structural blast analyses. The performance of both methods was also compared and highlighted by Sohn et al. (2013) based on pressure-impulse (P-I) diagrams.
Despite the technological advancements in FEA, the quality of outputs from these numerical analyses strongly depends on the skill and experience of the user. Appropriate modelling techniques, e.g. mesh densities, load and boundary conditions, material models, and element formulations, are essential to ensure that the FE model properly represents the real/actual structure. In structural dynamics, much has to be considered for the finite element (FE) formulations governing the parameters of interest for a particular engineering problem. Boh et al. (2004) recommended the use of first-order reduced integration shell elements for efficient and accurate FE simulations of blast response. Schwer et al. (2005) conducted a three-dimensional (3D) patch test based on the solid mesh proposed by Macneal and Harder (1985) to assess the performance of hourglass control functions via explicit finite element software, LS-DYNA. Sun (2006) demonstrated the compromise between reduced and full integration schemes for FE formulations in dealing with shear locking and hourglassing problems, whereby the details on both problems are as explained by Koh and Kikuchi (1987). As element formulations are intrinsically defined in commercial finite element software, their selection poses challenges that can directly influence the quality of the solution outputs.
The aim of the present study is to recommend relevant finite element (FE) formulations for blast simulation of a corrugated blast wall model by assessing the performance of the shell and solid elements with the aid of hourglass control functions in LS-DYNA. The outcomes of this study will provide greater acumen on the recommended FE formulations in terms of solution accuracy and computational cost, which can generally be applied in various FE software as most FE formulations are somewhat similar (Langer et al., 2017). Recently, Ng and Hwang (2017) conducted research on FE formulations on limited number of scenarios and extended results can be provided by the present study.

TYPES OF FINITE ELEMENTS IN LS-DYNA
In the present study, LS-DYNA is used to investigate the influence of relevant FE formulations on the structural behaviour of blast wall models subjected to explosive loading. Among various types of FE types, following three (3) FE types were only adopted in the present study.

• SECTION_SHELL • SECTION_SOLID • SECTION_TSHELL
The thin-shell (hereafter referred to as shell or used interchangeably with "shell" in this study), thick-shell, and solid finite elements have been preliminarily selected to model the corrugated panel and the supporting members of the blast wall.

Shell elements
As mentioned above, two types of shell element formulation such as thin-and thick-shell can be selected based on time efficiency, material type, simulation type (i.e., implicit or explicit) in LS-DYNA.  From previous studies, general features of abovementioned thin-shell elements were investigated by Stelzmann (2010) and Haufe et al. (2013). Generally, the thin-shell elements listed in Table 1 can be categorised as follows.
• Hughes-Liu shell formulation (EQ.1, 6, 7& 11) • Belytschko-Lin-Tsay shell formulation (EQ.2, 8 & 10) • Fully-integrated shell formulation (EQ.16) • Thickness enhanced shell formulation (EQ.25 & 26) Briefly, the Hughes-Liu shell formulation (EQ. 1) can be considered as a cost-effective solution as it is based on a degenerated continuum element formulation, in which 5-degree of freedom in local coordinate system and onepoint integration are adopted due to efficiency issues. It is also an effective method especially when large deformation needs to be taken into account. This formulation can also treat element warping. EQ.11 is also similar to EQ.1 except for the co-rotation system, such that EQ.11 requires additional computational cost. EQ.6 and 7 requires 3-4 times computational cost from adopting selective reduced integration (SRI) to avoid most hourglass modes.
The Belytschko-Lin-Tsay shell formulation (EQ. 2) is the default type in LS-DYNA, which is based on Reissner-Mindlin kinematic assumption (5DOF in local and 6DOF in global) and gives extremely cost-effective computational solutions. The bi-linear nodal interpolation with one-point integration is adopted. Fully-integrated shell formulation (EQ. 16), which is also based on Reissner-Mindlin kinematic assumption with 2×2 integration in the shell element plane, is recommended for implicit simulations. It does not degenerate to a triangle and requests 2 -3 times of additional computational cost, but with higher accuracy. The thickness enhanced shell formulation is also based on Reissner-Mindlin kinematic assumption with onepoint integration and bi-linear nodal interpolation (EQ. 25). In the case of EQ. 26, a 2×2 integration in the shell element plane is adopted while the Bathe-Dvorkin transverse shear correction helps to eliminate W-mode hourglassing. In addition, the linear strain through thickness feature is adopted.
The details of abovementioned four representative thin-shell element models are as described by Haufe et al. (2013).
The thick-shell elements shown in Table 2 can also be categorised as follows. This thick-shell element based on 8-node shell/solid is considered to be between thin-shell and solid element. The thin-thick shells (EQ. 1 and 2) are composed of 8-node shells with 2D stress state similar to that of thin shell. Basically, a penalty function is adopted to constrain the element thickness between top and bottom nodes. Once membrane stress is applied, only then can element thickness be changed. In general, the thin-thick shells depicted in Table 2 (EQ. 1 and 2) are not recommended due to efficiency issues in comparison to thin shells in Table 1.
Thick-thick shells (EQ. 3 and 5) also adopts the 8-node shell/solid, but is presumed to be under 3D stress state. In this case, the element thickness matter is resolved, which can be changed by thickness stress, however, it requires an extremely long computational time. In the case of EQ. 5, shear locking and hourglass issues are resolved and the laminated shell theory is applicable. This essentially helps to solve the engineering problem of bending with one element over thickness. It is developed for modelling thick composite structures whereby improper element ratio can also be considered. Table 3 shows the solid element formulations. There are several types of elements, however, only few are accentuated in the present study. Mesh-free hexahedron -A standard element (=EQ. 1) is set as the default, which consists of 8-node hexahedron solid element with trilinear shape functions. It adopts reduced integration, i.e., one-point integration in the middle of the element. A fully integrated element (=EQ. 2) is similar to the default element. This element adopts eight integration points which consumes 2-3 times additional computational cost than that of EQ. 1. It considers hourglass mode issues but may bring about shear locking and lower deformation problems.

Solid elements
A fully integrated quadratic 8-node element with nodal rotation (=EQ. 3) has 6-DOF with 14 integration points, which demands more expensive computational cost than EQ. 2. It is not listed in Table 3, but EQ. -1 and EQ. -2 can also be used, which were developed for solving the issues of fully integrated selective reduced hexahedron without shear locking. Normally, EQ. 1, 2, and 3 are recommended for modelling structures.

TARGET STRUCTURE
3.1 Blast wall design Figure 1: Configuration of the target blast wall model.
Blast wall structures are integrated installations on offshore topsides to minimise the effects of explosive loading. The present blast wall model was appropriated from a research report by HSE (2003), who provided relevant experiment data that are beneficial to this study. In addition, various studies on blast wall analysis and design optimisation have recently been conducted by several researchers (Kim et al., 2014;Hedayati et al., 2015;Syed et al., 2016;Li et al., 2017;Hao et al., 2017). The configuration and dimensions of the target blast wall consisting of a corrugated panel and connecting parts including angle, flexible angle, and I-beams are shown in Figure 1. In addition, this blast wall is a ¼ scale model of the real structure that was constructed and tested by HSE (2003). The material properties, i.e., material type, density, elastic modulus, yield strength, and Cowper-Symonds coefficients, are summarized in Figure 2. The Cowper-Symonds constitutive equation in Eq. (1) proposed by Cowper and Symonds (1957) is commonly applied in the ships and offshore industry (Park et al., 2015a;2015b;Choi et al., 2016;Kim et al., 2018) to consider the dynamic or strain rate effects, which are obtained via dynamic tensile tests. (1) where Yd  = dynamic yield strength, Y  = static yield strength,   = strain rate, Cand q = Cowper-Symonds coefficients obtained through curve-fitting. The dynamic characteristics of the materials are illustrated by the plot of dynamic yield strength normalised by static yield strength versus strain rate in Figure 2.

Applied blast loading
The applied blast loadings were idealised as triangular load curves by noting the rise time, r t duration time, d t and peak pressure, peak p as shown in Figure 3(a) as inputs into the numerical pre-processor. The maximum and permanent mid-span displacements are also demonstrated in Figure 3(b).

APPLIED EXAMPLE
In the present section, performance of pre-selected LS-DYNA element types, i.e. thin-shell, thick-shell and solid for modelling of the target blast wall were investigated through which relevant finite element (FE) formulations were selected for further numerical studies with recommendations.

Experimental test results
This section describes the experiment conducted by HSE (2003). Target blast wall were tested, by HSE (2003), in the pulse pressure loading rig developed at the University of Liverpool. Two (2) loading directions, i.e., positive "A" and negative "B" were considered for the pulse pressure tests, as shown in Figure. 4(a). In present study, totally seven (7) and three (3) pulse pressure load profiles were selected for two (2) opposing loading directions, namely positive "A" and negative "B", respectively, as clarified in Figure 4(b) and 4(c), respectively. These load profiles were used as load inputs for the numerical simulation while the previous blast test results, i.e. the maximum and permanent displacements summarised in Table 4. This data will directly be used for comparisons with the numerical simulation results in the following sections.

Numerical modelling
In ensuring robustness and safety of offshore structures, the possibilities of failure due to accidents such as explosion or fire should be anticipated and clearly reflected in the design stage. Due to high costs and time restraints, numerical methods such as nonlinear finite element method (NLFEM) are widely favoured in the offshore industry to ascertain structural responses, in contrast to experimental or destructive testing.
Taking advantage of symmetry, the FE model of the blast wall was simplified as one corrugation bay, half corrugation bay, and quarter corrugation bay. Figure 5 presents the displacement-time plots of these simplified FE models together with the full experimental model subjected to a peak pressure of 1.92 bar (loading scenario A7), to compare their permanent displacements with the test (HSE, 2003). Based on this investigation on the effect of model sizing, the quarter section model was adopted for subsequent analyses in this study, appreciating the reduced computational costs without compromising result accuracy. LS-DYNA explicit FE solver was used to perform the numerical simulation. The structural responses of the blast wall model subjected to a range of pulse pressure load profiles illustrated in Fig. 4 were evaluated. Pertaining to boundary conditions, the upper edge of the model was assumed to be fixed with both sides of the model set to be symmetrical in the transverse direction; the bottom edges were set to be symmetrical in the longitudinal direction, while a uniformly-distributed time-dependent idealised pulse pressure loading was applied all over the surface of the corrugated panel, as illustrated in Figure 6.
In the previous study by Sohn et al. (2012), 4mm of mesh size was adopted to the entire FE model. For the confirmation, mesh sensitivity analysis was conducted by adopting 2 mm, 4 mm, 8 mm, and 16 mm of mesh sizes in the present study shown in Figure 7. From the obtained results, we have confirmed that 4 mm of mesh size is relevant to be applied to FE modelling of the blast wall.  Material model No. 24 in LS-DYNA was used to represent the nonlinear dynamic behaviour of the structure by specifying the material strain rate parameters.

Assessment of FE types and FE formulations in LS-DYNA
This section is divided into two parts. First, the assessment and selection of FE types, i.e. solid, thin-shell, and thick-shell are addressed in section 4.3.1, from which further assessment of the selected FE types are discussed in the context of performance of FE formulations, i.e. reduced or full integration, and hourglass control in section 4.3.2. The structural responses, i.e. maximum and permanent displacements were validated against the test results (HSE, 2003). In addition, the obtained outcomes were discussed based on accuracy of numerical simulation results as well as computational cost.

Selection of FE types
Five (5) representative cases as shown in Table 5 were generated to study the blast response of target blast wall models based on combinations of pre-selected quadrilateral thin-and thick-shell, and hexahedral solid finite elements in LS-DYNA. Fig. 8 provides an overview of blast responses for extreme load scenarios in both loading directions, i.e. A7 ( peak p =1.92 bar) and B10 ( peak p =1.18 bar) for each of the five cases in terms of peak structural displacements, with respect to experimental measurements provided in Table 4.
Note: S, SH, and SHT denote solid, thin-shell, and thick-shell, respectively in present study.  (Table 5 can be referred to for the naming of each model).
Statistical approach was employed to analyse the FEA results. Fig. 9 shows the comparison between the present FEA and test results (HSE, 2003) for all loading scenarios and for all combinations of FE types with the calculated means, coefficients of variances (COVs), coefficients of determination (R 2 ) and standard error of the regression (S). Details on statistical analysis results can be referred to in Table A.1 to A.5 (Appendix A).  (Table 5 can be referred).
In Fig. 9(b), mean and COV values could not calculated because zero (0) deformation was measured from the testing. Therefore, FEA/test could not be done. In this regard, R 2 and S values were added on behalf of mean and COV values for the comparison. With regards to the accuracy of numerical simulation, models from Case 1 or Case 2 model may be suggested for further analyses based on the obtained outcomes as presented in Figs. 9(a) and (b). Furthermore, the performance of all cases can be sorted in the sequence of increasing computational costs as Case 1 (Cheap)< Case 2 < Case 3 < Case 4 < Case 5 (Expensive), while the details are referred to Table B.1 in the Appendix part. Throughout this study, the Intel® Core™ i7-6800K CPU @ 3.40GHz computer processor with 64-bit Operating system was used.
Therefore, in this section which covers "Selection of FE types", Case 1 and Case 2 can be recommended to users like structural designers for the analysis of offshore blast wall structure subjected to explosive loading. More details on Case 2 which apparently gives higher accuracy than Case 1 is scrutinised in the next section.

Selection of FE formulations
In section 4.3.1, shell and solid elements were recommended for modelling corrugated panel and connection parts, accordingly. In this section, the performance of several pre-selected FE formulations will be assessed. Table  6 shows four (4) FE models, comprising the combinations of reduced-and full-integration of solid and shell elements -systematically Types I, II, III, and IV -for detailed investigation. Each FE model is named according to the designation of FE formulations in LS-DYNA (refer to section 2) in the following sequence: supporting memberscorrugation panel-hourglass control.
For example, S1-SH16-8 refers to a model consisting of reduced integration solid supporting members and full integration shell corrugated panel with hourglass control EQ.8, while S3-SH2-0 refers to one consisting full integration solid supporting members and reduced integration shell corrugated panel without hourglass control. Table 5 is referred for the abbreviations. For the relatively thicker connection parts (supporting members), the number of through thickness solid elements can be another factor that influences the numerical results. Figure 10 demonstrates the relationship between several configurations of through thickness element distribution and the maximum midspan displacement (HSE, 2003), which clearly indicates that at least two layers of elements are required for a proper representation of the bending behaviour. Thus, for efficient FE simulations, configuration (c) with 3, 2, and 2 layers of elements through the thicknesses of I-beam, flexible angle, and angle, respectively, has been adopted for all the following FE models in this study. Through comparisons with the test results (HSE, 2003), it was observed that the influence of FE formulations rises exponentially with increasing peak pressure. Figure 11(a) compares the displacement-time histories of all models (Types I to IV) subjected to 1.92 bar peak pressure (loading scenarioA7), which clearly indicates the capabilities of these FE formulations in predicting the dynamic responses of the blast wall model.
The excessively overestimated responses from Type I (S1-SH2-0) and Type II (S1-SH16-0) models were due to the hourglass modes that generated zero energy in the affected solid elements, particularly at regions of large deformation, causing the connection angles to lose stiffness hence exaggerating the maximum response. The hourglassing phenomenon is shown in Figure 11(b). While reduced integration (RI) hexahedral solid and quadrilateral shell elements are prone to hourglassing, Type III model (S3-SH2-0) was effective in mitigating the undesirable elemental "defects" in most cases that were subjected to low peak pressures.
For instance, the permanent displacements in loading scenario A7 was slightly underestimated by the fullyintegrated (FI) solid elements, which might be due to the shear-locking phenomenon that over-stiffens the responses. In addition, the effect of number of integration point, which caused the different outcomes, should also be carefully taken into consideration; the hourglass mode can occur if only one (1) integration point is adopted for the FE simulation. Furthermore, relevant additional options may be required to prevent hourglass mode.
In contrast, Type IV (S3-SH16-0) model, which consists mainly of FI elements, has shown little deviation from that of Type III model, implying the insignificant influence of shell element formulations on solution output. Models without hourglass control at the connection (solid elements) and the corrugated panel (shell elements) during maximum response are shown in Figure 11(b). Since hourglassing was observed only in the RI solid elements of S1-SH2-0 and S1-SH16-0, a simple deduction can be made such that RI solids would suffer the numerical shortcomings, in which they tend to be excessively flexible, hence the overestimation of permanent displacement. However, hourglassing did not occur in the FI solids shown in Figure 11(b) -connection part modelled by solid element. Figure 11: FEA outcomes of Type I-IV models without hourglass control (loading scenario A7).
Essentially, hourglass control functions can be viewed as correction terms for numerical integrations, which introduce internal nodal (hourglass) forces that are proportional to the components of nodal velocity or nodal displacements to counterbalance the zero-energy modes (LS-DYNA, 2014). The available hourglass control functions in LS-DYNA were applied in accordance to Table 6, to all models. Basically, under-integrated solid and shell elements are prone to hourglass modes, which can be handled in two ways, e.g. by applying hourglass control functions in LS-DYNA and/or by mesh refinement. Figures 12(a) and (b) present the extent of hourglassing in terms of hourglass energies generated in the solid and shell elements, respectively, with respect to element mesh size of a representative analysis case. It is evident that the hourglass energy is directly related the fineness of element mesh, though this aspect of modelling is often insignificant for shell elements due to the location of integration points. As a rule of thumb, the generated hourglass energy in an element shall be well below one-tenth (or 10%) of the total energy generated in that element. Findings in Figure 12 agreed well with the adopted mesh size as shown in Figure 7. From Figures 13(a) and (b), Type I and Type II models both have seen result improvements for the RI solids with the aid of hourglass control functions, whereas Type III model in Figure 13(c) was not improved by this addition of corrective forces as FI were already in use. Based on the overall results shown in Figure C.1 (Appendix C), including all ten (10) load scenarios (A1 to B10), Type III and Type IV models are not satisfactory in predicting the structural response at high peak pressures as they tend to be excessively stiff, thus underestimating the permanent displacements. Figure D.1 in Appendix D presents the effective plastic stress and strain distribution contours for all Type I-IV models subjected to loading scenario A7. Due to advantages in computational efficiency and ability to overcome shear locking, RI schemes are widely implemented in explicit FEA codes. However, the downside of these elements is their tendency to introduce hourglass modes of deformation, in which neither stresses nor strains are generated in the affected elements, thus illdefining the resulting structural response.
While FI elements are effective in dealing with hourglass instabilities, their major drawback is, as opposed to that of RI elements, over-stiffening of the responses by shear-locking. In short, shear locking and hourglassing are two compromised factors between the two integration schemes. Thus, as it may be difficult to eliminate hourglassing, some form of hourglass control is required in the FEA.
As in the previous section, statistical analysis issued to compare the present FEA results against test results (HSE, 2003) for both maximum and permanent displacements for Case 2 only, i.e. four combinations of solid and shell elements according to Table 6 As S1-SH2-0 and S1-SH16-0 were affected by hourglassing, their calculations are unaccredited and the values shown in Figures 14(a) and (b) were discarded. From Figure 14(a), Types I and II with slightly higher R2-values performed well in predicting maximum displacements compared to slightly underperformed Types III and IV, while from Figure 14(b), Types III and IV have somewhat outperformed Types I and II in predicting permanent displacements. However, though with decent capabilities, Types III and IV models were deemed too conservative in dealing with cases involving high peak pressures, in addition to requiring much higher computational costs as shown in Table B.2 (Appendix B).
Type IV (Expensive) < Type III < Type II < Type I (Cheap) Comparing the computational cost, the same PC has been used as mentioned in Section 4.3.1 and the following order is defined. Thus, the use of RI elements with appropriate hourglass control function is indeed recommended for cost-effective solutions.

CONCLUDING REMARKS
In the present study, the influences of solid and shell FE formulations with additional hourglass control functions on the maximum responses of corrugated blast wall has been investigated using LS-DYNA explicit finite element solver. It is wise to take advantage of thickness variation over the entire target structure when selecting the representative finite element (FE) types, i.e. thin-shell, thick-shell or solid to represent the model parts.
The selection of element integration scheme, similar to other considerations in FE modelling, is essential to the quality of FEA solutions. Reduced integration (RI) elements are favourable in explicit dynamics analyses, given its high speed and robustness under high structural distortions, whereas full integration (FI) elements are more typical in implicit analyses. Although FI solid elements perform consistently well in predicting maximum responses under low peak pressure, they are very costly compared to their RI counterparts.
Furthermore, FI solid elements are not suitable for predicting responses of high peak pressure as the effect of shear locking increasingly falsifies or over-stiffens the responses, which had been observed through the comparison with experimental measurements. Moreover, the number of integration points should also be carefully taken into consideration in order to prevent the shear locking phenomenon of solid element.
In contrast, RI solid elements associated with relevant hourglass control functions can be used to obtain satisfactory estimations of blast responses for all peak pressures in much shorter computation times. The present recommendations are presumably software independent and generally govern different FE software.
In the future, researchers may want to delve into the effect of shear locking phenomenon as it should be further studied by including solid element formulations -1 and -2 in LS-DYNA.