Comparative study of shell element formulations as NLFE parameters to forecast structural crashworthiness

1 Introduction

The potential for massive casualties involved following maritime accidents, e.g., collision, grounding, and explosion, means that impact loading must be considered to avoid fatalities [1–5]. Estimation of the structural performance when subjected to this loading type has developed since the first damage energy formula was introduced by Minorsky [6], and refined versions of this formula have been presented by Woisin [7], Vaughan [8], and Zhang et al. [9]. Various studies have followed this finding, i.e., through the improved estimation obtained by including data from analytical analyses [10,11] and laboratory experiments [12,13]. On the other hand, the rapid development of computational technology that took place in the late 1990s provided the finite element (FE) approach as one alternative calculation method for structural investigation, especially for modeling delicate subjects and characteristics. An advantage of this methodology is the reduced risk of time loss due to rejection or failure, which is the main reason that engineering phenomena are suitable to be idealized using the FE approach. Nevertheless, considering the complexity of the nonlinear behavior of materials and time dependency during impact, a detailed definition of the material and other parameters in terms of the geometry and boundary conditions needs to be developed using the nonlinear finite element analysis (NLFEA) [14–17].

Designs of impact loading for application in maritime accident investigation are created by considering physical parameters, e.g., impact speed, interaction angle, vessel mass, and target location, which involve model idealization and assumption, i.e., calculation parameters in NLFEA, which inevitably lead to diverging results; proper quantification is required to control the validity of the method. One of the fundamental parameters is element formulation (EF), which has to be defined in every FE calculation. This formulation is embedded in the geometrical model (marine vessels are mostly assumed to be thin-walled structures [18–20]), which directly contributes to the magnitudes and contours of structural performance. Although its importance influences the result estimation, the formulation has not been widely discussed as studies related to this subject are scant.

In the realm of engineering, finite element analysis (FEA) has proven to be an invaluable tool for simulating and understanding the complex behaviors of materials and structures. This methodology allows engineers to assess the performance of various systems under different conditions, aiding in the design, analysis, and optimization of structures subjected to diverse loads. However, the nonlinear behavior of materials and the temporal aspect of progressive phenomena, especially during impact scenarios, pose significant challenges that demand a refined approach. In this review, we delve into the advantages and challenges of using NLFEA for investigating maritime accidents, with an emphasis on the impact loading dynamics.

This study aimed to undertake a comparison of a series of impact loading scenarios, each grounded in a geometrical assumption. These scenarios will be calculated using NLFEA in order to assess the influence of EF on the resulting outcomes. The initial phase involves benchmarking against pioneer laboratory experiments and an established shell formulation. The primary objectives here are to validate the reliability of the current NLFEA methodology and make preliminary predictions regarding the effects of different formulation types on the structural performance. Subsequently, an expanded investigation will be conducted. This phase will involve a modeled double-hull passenger vessel subjected to ship–ship collision conditions. The assessment will consider performance criteria, such as energy, effective strain, and displacement, to precisely quantify how the chosen shell EF impacts the structural response of the ship under impact loading conditions.

2 Shell formulation and theory

The application of a shell element in idealizing thin-walled structures using the FE approach has conventionally occurred due to hardware limitations and efficiency of analysis. Thin-walled structures in the real world are one type of three-dimensional shell body, which can be modeled via three different methodologies, as presented in Table 1. Implementation of the shell method, a component of shell formulation, requires good understanding. The evolution of shell models began with the Kirchhoff–Love theory regarding membrane and bending effects, which was introduced in 1888; two points are underlined as the main topic, i.e., membrane and bending action and inextensible deformations [21]. The fundamental assumption of this theory is that the cross sections remain straight, unstretched, and normal to mid-surface. Contradiction occurred and required a modification of the material law when stress σ _zz and ε _zz remain zero. The next idea, proposed in the middle of the 20th century, the Reissner–Mindlin–Naghdi theory, proposed transverse shear strain. This theory excluded a criterion of the previous theory, i.e., cross sections normal to mid-surface [22]. Nevertheless, the weaknesses of the Kirchhhoff–Love theory had not been solved. The Ahmad–Iron–Zienkiewicz theory from 1968 introduced a concept, the so-called degenerated solid approach, in which shell theory is defined as semidiscretization of a 3D continuum [23]. After 1990, the widely used assumption was 3D-shell FEs, solid shells, and surface-oriented (continuum shell) elements [24–26].

Table 1

Thin-walled structure characteristics

Original form	Method	Definition process
3D shell body (continuous)	Shell	Reduction into 2D (semi-discretization e.g., linear in ζ)	Shell model (discrete in ζ, continuous in ξ η)	FE discretization 2D shape functions (e.g., bilinear in ξ η)	6 degrees of freedom (DOF) per node, e.g., 4 nodes per element
	Solid	FE discretization 3D shape functions (e.g., trilinear in ξ η ζ)			3 DOF per node, e.g., 8 nodes per element
	Solid shell	Special EF, e.g., assumed strain formulation

The derivation of these theories’ evolution is implemented in FE calculation in the form of shell EF, which must be defined. Several formulations embedded in ANSYS LS-DYNA are described in Table 2, including Belytschko–Tsay (BT) (default), Belytschko–Wong–Chiang (BWC), Belytschko–Leviathan (BL), and fully integrated (FI) BT [28–31].

In recent years, there has been growing interest in the analysis of shell structures based on higher-order theories that can account for the dynamic effects and layer-wise (LW) variation of the field variables [32–35]. Among these theories are the LW approach, the equivalent single-layer (ESL) approach, and the equivalent layer-wise (ELW) approach. Some examples of these theories are the generalized higher-order LW theory for the dynamic study of anisotropic doubly curved shells with a mapped geometry; an ELW approach for the free vibration analysis of thick and thin laminated and sandwich shells; and LW and ESL models for smart multilayered plates.

3 FE benchmarking

Benchmarking in FEA is the process of comparing the performance and accuracy of different numerical models of a physical system. A benchmark FE model is a reliable model that can fully and correctly reflect the real behavior of the structure and has been verified by field tests. It is also important for validating the assumptions and approximations made in the finite EF, such as the material properties, the boundary conditions, and the mesh size, as well as for identifying the strengths and weaknesses of different methods, such as linear or nonlinear, static or dynamic, and implicit or explicit. Benchmarking can also help to establish best practices and standards for FEA in various fields related to engineering and applied mathematics.

3.2 Benchmarking results

The obtained results for crushing force in Figure 1 indicate that the current methodology produced good similarity. Both results displayed increment tendency with a maximum state recorded in the range of 175–200 mm; after this point, the force was reduced as a fracture occurred on the target plate. The setting of material characteristics, i.e., fracture/failure criteria, is also compared in this discussion, and a recent work by Cerik et al. [37], who proposed the Hosford–Coulomb criterion for shell structures, is considered.

Figure 1

Comparison of the test by Alsos and Amdahl [36] with the pioneering work in FEA by Cerik et al. [37].

The comparison shows the correspondence of the current work with the original test and a pioneering FE analysis in which the maximum state was reached at 175–200 mm, but the results of the study by Cerik et al. [37] showed early analysis termination as no response was recorded at a displacement of 225 mm. A variation in mesh size is considered a major influence in FE analysis and requires calibration to obtain satisfying results. The tendency in Figure 2 illustrates that in terms of fracture occurrence (indicated by the maximum state of force), the mesh size of 15 mm in the current work is the closest to the results using the Hosford–Coulomb criterion [38–40] with a displacement of 200 mm. There is also similarity of these results in terms of the fracture pattern, as shown in Figure 3, which reveals that crescent-moon-like tearing took place during the impact on the plate panel. Based on benchmarking using laboratory tests [36] to validate the reliability of reproducing the actual test/phenomenon, and FE simulation [37] to verify material fracture characteristic in this stage, the satisfactory nature of the current FE methodology was confirmed.

Figure 2

Influence of mesh size on crushing force, accounting for the variety of shell formulations.

Figure 3

Structural contours of the test that compared the pioneering work and the current work to calibrate fracture behaviors: (a) laboratory experiment using Alsos and Amdahl [36], (b) pioneer FE using Hosford-Coulomb Cerik et al. [37], and (c) current FE using laboratory experiment setting.

4 Methodology and analysis preparation

An extended study is performed in the following section to assess the responses of larger structures, i.e., ships. Ship–ship collision is assumed as the main scenario and involves a double-hull passenger ship as the target ship and a cargo ship as the striking ship. Regardless, a previous study has shown that ship–ship collision can be prevented using radar, but the fact is that these severe accidents still happen all over the world [41]. The target ship is idealized in FE analysis to have deformable thin-walled structures [42–45], with three regions set to apply proper mesh sizes. The smallest size, 8–9 mm, was applied to the core region, which was expected to experience the direct impact of the striking ship, with major deformation expected to occur. The 9–10 mm element size was embedded on the transition region as a hub between the core and the outermost regions. A larger size with size 10 mm and above was implemented in the outermost zone as it is predicted that no fracture will occur or expand into this region. The striking ship, on the other hand, is designed as a rigid structure, which means that all the absorbed/internal energy in the collision is focused on the target ship. The striking ship is given a uniform velocity, approaching the targets at 15 kts or approximately 7.72 m/s. As impact occurs between the two ships, the contact model must be defined with the impact that was modeled by automatic contact and the friction coefficient between steels implemented.

The conditions of the impact were dominated by the movement of the striking ship, while the target ship was restrained for its axial and rotational displacement at its original location. The fixation was applied at the ends of the deformable structure, including web frames and bottom members. The schematic illustration presented in Figure 4 indicates the three vertical locations on the target ship selected as impact points for the striking ship. Fundamentally, these points were determined based on the possible increase in the striking ship’s draught. Precisely, these were located on the middle deck (approx. 9 m from the baseline (BL)), side shell (approx. 10 m from the BL), and main deck (approx. 11 m from the BL) and labeled Points 1–3, respectively. In order to quantify the effect of EF on crashworthiness criteria and time processing, the structures were applied by four augmented shell EFs: the BT EF = 0, the BL (EF = 8), the BWC (EF = 10), and the FI version of the BT (EF = 12) [28,29].

Figure 4

Designed collision scenario for the extended study of formulation influence on crashworthy structures.

5 Extended study of structural crashworthiness

The physical configuration in scenario design refers to the arrangement of the target locations that are subjected to impact loading. The target locations can be chosen based on the expected collision scenarios or the critical regions of the ship structure. The results of the combined parameters, such as EF, material model, and physical configuration, are assessed to quantify structural performance under impact loading. Structural performance can be measured by various crashworthiness criteria, such as energy absorption, which is the amount of kinetic energy lost by the structure during deformation; effective strain, which is a measure of the average plastic deformation in a region of interest; displacement, which is the change in position or shape of the structure due to loading; and damage extent, which is the size and location of cracks or fractures in the structure. These criteria can indicate the level of deformation, fracture, and failure of the ship structure due to impact loading.

5.1 Internal energy and crushing force

During impact, the involved structures produce a response in the form of energy as a consequence of the destroyed or deformed members. This energy type is called internal energy or absorbed energy in the literature related to crashworthiness [46–48]. The results displayed in Table 3 indicate that quite a different energy level is shown by the BL formulation, with an estimate of 13% compared to the highest energy produced by the BT formulation. Kinetic energy is another energy response that expresses the movement of the involved parameters. The overall tendency to give a prediction of higher internal energy is caused by higher kinetic energy. A change in the hourglass energy reflects irregularities during the structure’s deformation and ultimate destruction. Lower hourglass energy is desired in the analysis, especially in NLFEA. The FI version of the BT is proven to be the most powerful formulation that is capable of fully suppressing this energy. The total energy is the sum of all energy types in the collision process. An assessment is conducted on internal energy, kinetic energy, hourglass energy, and total energy and finds that other energy, in the range of 14–17%, also contributes to the process. This energy is defined in the FE analysis as contact (sliding) energy, rigid wall energy, and system damping energy. Overall assessment in terms of energy shows good similarity between all the selected formulations, of which the BL is noted as the formulation that underestimated internal energy compared to the others. On the other hand, the reliability of NLFEA was successfully guaranteed by the FI version of BT, followed by BL with the lowest value for the hourglass energy/internal energy ratio.

Table 3

Energy results of the selected shell formulations

Formulation	Kinetic	Internal	Hourglass	Total	I/K Ratio	H/I Ratio	I + K + H	(I + K + H)/T	(I + K)/T	Other energy
	K	I	H	T
	(J)	(J)	(J)	(J)	(–)	(–)	(–)	(–)	(–)	(–)
BL	1,492,730	5,982,160	148,742	9,233,940	0.2495	0.0996	7,623,632	0.8256	0.8095	0.1744
BT	2,204,380	6,875,120	493,751	12,105,200	0.3206	0.2240	9,573,251	0.7908	0.7500	0.2092
BT-FI	1,425,950	6,787,900	0	10,148,000	0.2101	0	8,213,850	0.8094	0.8094	0.1906
BWC	1,475,710	6,815,060	389,818	10,210,200	0.2165	0.2642	8,680,588	0.8502	0.8120	0.1498

The crushing force is the acting force at the moment the deformable structure experiences a fracture and the latter is destroyed. The tendency in Figure 5 is very similar until bow displacement at 0.25 m, when it can be predicted that a major fracture of the target ship has not been manifested. After this displacement distance, the BL formulation shows a higher response, while for the BT, the FI version of BT, and the BWC, it takes place approximately until the end of the collision process. After displacement at 1.25 m, the force of the BL was below that of the other shell formulations. Nevertheless, considering typical fluctuations and patterns, the results can be considered identical, with no significant discrepancy.

Figure 5

Tendency of the crushing force during ship–ship collision.

5.2 Strain and stress contours

Strain is indicated in nonlinear FEA when material and structural elements experience deformation, fracture, and finally destruction. The behavior of strain presented in Figure 6 represents an element on the upper side of the tearing length (outer shell) located in the impact region. The results indicate that the element had not fully experienced fracture as the ultimate strain was set to be 0.2. The BWC formulation produced final strain in the range of 0.075–0.1, which was also shared with the strain contour in the inner shell, as displayed in Figure 7. Interestingly, the mentioned formulation and the BT showed a longer strain contour, which indicates destroyed structures due to initial penetration on the outer shell lending more significance to the inner shell compared to the BL and the FI version of BT.

Figure 6

Behavior of failure strain on the side shell during impact.

Figure 7

Failure strain of the selected shell EFs.

The similarity of the tendency of stress behavior (Figure 8) shows that the initial impact did not produce any fracture until displacement at 0.25 m. Beyond this point, shell formulations tended to fluctuate in the same pattern until 1 m. The stress level experienced at the end of the calculation was highest with the BL formulation, while other formulations were reduced below 2 × 10⁸Pa. This phenomenon possibly occurs since the selected element may interact with other impacted members, and the latter causes higher fluctuation in stress history, while a lower stress level indicates that beyond a certain bow displacement, the selected element on the side structures of the target ship was unaffected (or affected only at a low level) by the collision process. The comparison is more explicitly assessed at a global scale as presented in Figure 9, in which the longer strain contours of the BT and BWC formulations also appear as a longer stress contour with quite a high level on the inner shell, especially near the bow structure of the target ship.

Figure 8

History of critical stress on the selected spot above the tearing damage on the side structure.

Figure 9

Critical stress contour of the target ship under ship collision.

5.3 Accumulated displacement

Displacement in this discussion refers to the component/member being displaced from its original position due to ship collision. The results in Table 4 are arranged by the displacement direction for each member, with the most displaced member during impact being the middle deck in the transversal direction at a range of 0.1577–0.3908 m. The car deck was found to be the least affected by impact as all formulations produced very similar results in terms of displacement to the longitudinal, transversal, and vertical directions. Nonlinearity was hardly observed for this member, which indicates that the crushing of other members did not significantly affect the position of this member. The overall results indicate that the displacement level of members produced by the BT, the FI version of BT, and the BWC formulations was the most similar; however, in the case of the structural frames, the BT and the BWC were very similar in terms of displacement in three directions. Nevertheless, quantifying the formulation effect in specific ways, e.g., on member behavior, requires extensive observation due to the nonlinearity and dynamic characteristic of structures under impact loading.

Table 4

Displacement of the involved member in ship collision

Direction	Member	Formulation
		BL	BT	BT-FI	BWC
Longitudinal (m)	Middle deck	0.0484	0.1041	0.0631	0.1129
	Car deck	0.0003	0.0003	0.0003	0.0003
	Side shells	0.0251	0.0455	0.0351	0.0494
	Struct. frames	0.0193	0.0408	0.0269	0.0463
	Upper/main deck	0.0006	0.0011	0.0007	0.0016
Transversal (m)	Middle deck	0.1577	0.3510	0.3908	0.3908
	Car deck	0.0001	0.0002	0.0002	0.0002
	Side shells	0.0997	0.1629	0.1789	0.1789
	Struct. frames	0.0826	0.1521	0.1692	0.1692
	Upper/main deck	0.0005	0.0018	0.0027	0.0027
Vertical (m)	Middle deck	0.0569	0.0621	0.0220	0.0544
	Car deck	0.0075	0.0077	0.0058	0.0076
	Side shells	0.0107	0.0128	0.0107	0.0128
	Struct. frames	0.0052	0.0158	0.0051	0.0155
	Upper/main deck	0.0075	0.0218	0.0061	0.0380

As a solution, a larger-scale observation of the structural displacement contour was conducted as presented in Figure 10, which explicitly documents the behavior of two important members of the side structure, i.e., the outer and inner shell. There was an indication of similarity in member behavior in terms of the overall contour for the BT and the BWC formulations, whereas a wide range of displacement on both shells was found from their produced results. The FI version of the BT actually shared an identical pattern with the BL, but the displacement on the inner shell of the latter formulation was found to be very small. This identification process regarding the displacement behavior provides a good correlation with the energy response, for which the BL produced the lowest internal energy and total energy.

Figure 10

Displacement level of the side structure involved in ship collision.

5.4 Hourglass behavior and time analysis

Hourglass modes are nonphysical, zero-energy modes of deformation that produce zero strain and no stress. Hourglass modes occur only in under-integrated (single integration point) solid, shell, and thick shell elements. This energy is undesired by researchers, especially in the NLFEA field, as it reduces the accuracy of calculations. To confirm that the nonphysical hourglass energy is small, the energy should be relatively small in comparison with the peak internal energy for each part/member – less than 10% as a rule of thumb [46].

An initial assessment of the global results from the previous section (referring to Table 3) leads us to conclude that only BL and FI BT produce very good results in terms of the hourglass-to-internal-energy ratio, with values of 9.96 and 0%, consecutively. However, the hourglass assessment addresses each part of the structure, which produces calculation results for each member directly involved in impact loading. Figures 11–15 record the hourglass energy of the target ship during ship collision. The results suggest that all shell formulations are acceptable for NLFEA. This is also clear from the analysis of ships during collision, performed in previous studies and examined in considerable detail [49–54].

Figure 11

Hourglass pattern during collision of the local member: car deck.

Figure 12

Hourglass pattern during collision of the local member: structural frames.

Figure 13

Hourglass pattern during collision of the local member: main deck.

Figure 14

Hourglass pattern during collision of the local member: middle deck.

Figure 15

Hourglass pattern during collision of the local member: side shells.

Another important parameter in FEA is time simulation, as described by Bathe [55]. There are also several considerations when performing the FEA: for instance, the mesh size selection [56–58] and failure criteria [59–62]. The common material properties of steels used in FEA models can be found in the study conducted by Ridwan et al. [63]. Reliability and effectiveness are two different but related aspects to consider in FEA; reliability is an extension of accuracy, while effectiveness represents an acceptable time frame. Deploying a very fine mesh or FI version formulation is a way to ensure that this parameter is fulfilled. However, the time of a simulation can increase, and, in some cases, leads to an unacceptable processing time. On the other hand, a very coarse mesh or selection formulation with a rapid calculation ability may indeed have good effectiveness, but the accuracy of the results could be questionable [64–70]. Based on these formulations, processing time must be assessed as a follow-up of the previous discussion focused on the accuracy of results (reliability). Table 5 suggests that the fastest results were produced by the BL and BWC formulations. Two versions of (BT and BT-FI) produced the slowest results. Although the BL formulation has a good hourglass-to-internal-energy ratio, for an analysis involving more complex structures and arrangements, the FI-BT is still recommended as it is capable to ensure hourglass to be suppressed.

Table 5

Time processing of the collision analysis

Information	BL	BT	BT-FI	BWC
CPU time per zone cycle (ns)	156	116	277	147
Problem cycle (cycle)	325,226	778,499	323,327	42,736
Total CPU time (s)	3,767	6,719	6,637	4,463

6 Conclusions

This work was addressed to assess the effect of shell EFs on structural crashworthiness in NLFEA. Calibration of the current FE methodology was conducted in the initial stage of this research by considering a pioneer impact test on a stiffened panel. The results indicated good similarity between the current results and the laboratory test. This was confirmed by a previous work that used FEA to simulate the same impact case to calibrate fracture behavior. Satisfactory results were achieved in this stage, with a mesh size of 15 mm of the current configuration producing the most similar results to the previous simulation.

Extended study was performed on a more complex structure, i.e., a double-hull structure, with the objective of quantifying the effect of shell formulation on the structural response produced. With regard to energy, hourglass was proven to be entirely avoided by deploying the FI version of BT, while the ratio of global energy indicated that the BL was also a viable formulation for NLFEA. A close correlation between the internal energy and stress–strain produced lower energy and caused a lesser stress level and smaller strain contours. Crushing force, on the other hand, showed homogeneity across all formulations, producing a result when major damage on the side structure had not occurred or was in the range of approximately 0–1 m. Based on the history of critical stress estimated by the two versions of BT and BWC, in the end, they produced stress reduction compared to the middle of the impact process, which indicated that the selected element was no longer mainly affected by the crushing process. For evaluating the hourglass energy, it was important to ensure a ratio of hourglass-to-internal energy per part of 10%. All formulations matched the criterion and were capable of being deployed as part of the structural configuration in NLFEA. Nevertheless, the FI-BT is highly recommended for use as a shell formulation since the ratio in the global scale had an hourglass value of zero. Moreover, in terms of reliability, time processing as part of effectiveness had to be assessed as a consideration in the NLFEA process. The BL formulation produced the fastest results, followed by the BWC. On the other hand, the FI BT required more processing time as compensation for the suppression of hourglass energy.

Calibration and benchmarking are important aspects in FEA, especially for analyzing highly dynamic and nonlinear phenomena. Quantification of the formulation effect on the structural response at a very specific level, i.e., structure members, is exhausting work, so a rapid estimation is recommended to consider the global response. As a future extension of this study, one might consider the effect of shell formulations in calculating major and minor collisions.