Comparison of numerical approaches for structural response analysis of passenger ships in collisions and groundings

The dynamic response of ships following grounding and collision accidents may be influenced by structural topology as well as operational and environmental conditions. Traditionally, the consequences of such events may be assessed by crude empirical methods or laborious experiments. Computational methods offer a useful alternative in terms of accurately capturing crushing mechanisms also accounting the influence of surrounding water. This paper presents a benchmark study that compares the structural dynamic response by explicit nonlinear FEA approaches and the semi-numerical super-element method. Simulations for typical accident scenarios involving passenger ships confirm that implementing the influence of hydrodynamic restoring forces in way of contact may be useful for either collision or grounding. Yet, for grounding scenarios, the damaged area resulting from analytical simulations appears to be sensitive to the failure strain values adopted to model the rupture of the ship bottom floors.


Introduction
In recent years the continuous increase in international seaborne trade lead to highly intensive traffic activities that inevitably increase the likelihood of ship accidents (see Fig. 1). According to the European Maritime Safety Agency (EMSA), cargo ships are the most vulnerable ship segment (44.0%), followed by passenger ships (24.8%) [1]. For the period between 2014 and 2019 almost half (approx. 43%) of accidents related with collision, grounding and contact events. Historically, such events result in significant losses of human life as well as cargo and environmental damages (e.g. Ref. [2]. To mitigate the impact of consequences it is essential to design against Accidental Limit States (ALS) [3].
Fundamental to this is the development and validation of numerical and analytical methods. For example, Hong and Amdahl [4]; Lui et al. [5] and Zeng et al. [6] investigated the structural response of plated girder, plate or partial stiffened panel subject to the impact load by experimental tests. Simonsen [7] and Zeng et al. [6] developed an analytical method to predict grounding damage.
Pedersen and Zhang [8,9]; Paik and Seo [10] and Hong and Amdahl [4] introduced empirical formulae to calculate damages and structural crashworthiness following collision and grounding events. Glykas and Das [11]; Zhang and Suzuki [12]; Haris and Amdahl [13] and Nauyen et al. [14] presented Finite Element (FE) modelling procedures and carried out sensitivity analyses. More recently, Ehlers et al. [15] validated ship collision mechanics by large-scaled experimental investigations accounting for the influence of the added mass. Similar studies by Le Sourne et al. [16,17]; Yu and Amdahl [18]; Rudan et al. [19] demonstrated the importance of restoring forces on dynamic response during collision events. In an attempt to unify crashworthiness modelling assumptions, Ringsberg et al. [20] presented a benchmark of nonlinear FEA against experimental tests for side collision and grounding with heave velocity. In a similar fashion Brubak et al. [21] conducted a benchmark study on the influence of key FEA modelling parameters for the case of grounding against Simonsen [7]. Recently, Kim et al. [22] demonstrated the significance of the combined influence of evasive actions and environmental conditions (e.g. effects of hydrodynamic restoring forces) in way of collision and grounding. Paik et al. [23] and Youssef and Paik [24] provided statistical information on ship-ship collision and grounding with forward speed. This is useful information especially for structural safety assessment.
In principle, the validation of collision and grounding techniques demand comparisons against experiments that suitably idealise interactions between external and internal mechanics. Experimental studies are presented by Zhu et al. [25] and Zhou et al. [26] for grounding and more recently by Xu et al. [27] and Wang et al. [28] for collision. This paper presents a computational benchmark study that highlights the importance of realising uncertainties of traditional FE and contemporary super element models accounting for the impact of surrounding water on ship collisions and groundings. The methodologies presented are considered significant within the context of ship design for safety, i.e. the implementation of crashworthiness methods in future ship stability standards [29]. Tables 1(a) and 1(b) summarise the participants of this benchmark and their role.

Target structures
Two passenger ships namely Ship A (corresponding to a passenger ship) and Ship B (corresponding to a RoPax) were used in this benchmark (see Fig. 2). Table 2 presents their general particulars. Ship A has been used for collision simulations (as the struck ship) and grounding simulations. Ship B was the striking ship utilised in the collision benchmark only.

Shipship collision scenarios
In the collision benchmark, four scenarios were considered (see Table 3). The collision speed of reference scenario 1 originally implemented by Kim et al. [22] is based on the collision risk distributions introduced by Paik et al. [23]. In this scenario struck Ship A is assumed stationary. A detailed definition of parameters is shown in Fig. 3.

Grounding scenarios
For ship grounding, 5 different scenarios were considered. Those depend on the shape of rock, grounding location and grounding speed (see Table 4). The reference grounding scenario was selected based on the historical accidental data presented by Ref. [24]. Fig. 4 shows the definition of rock shape and grounding locations. According to SOLAS [31]; 2 m of height for both rocks may be considered acceptable. The paraboloid rock was modelled as follows: where C and E are coefficients (C = 3.7, E = 6). With 2 m of rock height, the longitudinal radius is 0.735 m and transverse radius is   0.577 m.

Effects of operational and environmental conditions
As demonstrated by Kim et al. [22] the effects of surrounding water (fresh or sea water) and associated environmental conditions may influence the impact velocity in way of contact, the resistance force and hydrodynamic actions (restoring force, added mass and damping force). Table 5 summarises these effects. Table 6 shows a comparison of the applied operational and environmental conditions in collision simulations.
The AALTO model accounted for the collision velocity following evasive manouvre as per Taimuri et al. [30] (see Table 7). The manourvring code by Taimuri et al. [30] uses a 4th order Runge-Kutta numerical integration to account for the ship equation of motions at each time step. STAR CCM+ [32] Reynolds-Averaged-Navier-Stokes (RANS) solver was used to calculate velocity-dependent resistance forces for the struck ship in sway and the striking ship in surge (see Fig. 5(a)). Fig. 5 demonstrates an example of the velocity-resistance force by STAR CCM+, pre-calculated time-velocity and applied time-resistance forces histories in FE analysis for scenario 2. In this case, the time-resistance force histories applied on the ship structure ( Fig. 5(c)) were converted to time-velocity curves ( Fig. 5(b)). It is noted that each scenario had different time-velocity/resistance force histories depending on the velocity during events. The detailed effects of initial impact velocity in 6 DoF, resistance force and hydrodynamic properties are explicitly explained in Kim et al. [22].
In AALTO, BV/ICAM and MSRC methods, hydrodynamic boundary conditions have been respectively applied on FE and superelement models via MCOL. This solver accounts for i) mass matrix, ii) hydrodynamic restoring forces, iii) water added mass, iv) buoyancy parameters and v) wave damping forces as per [17] and Hydrostar [33]. Accordingly, seakeeping parameters were calculated in 20 frequency steps (0.1-2.0 rad/s) in infinite water depth assuming zero forward speed for the struck ship and a forward speed of 5 or 10 knots for the striking ship.
MARIN used MarcolXMF [34]. This solver is based on a 6 DOF time domain rigid body hydrodynamics solver [35]. Restoring forces were computed by pressure integration over the hull geometry. Added mass and damping contributions were based [36]. This solver accounts for forward speed effects as well as the influence of hydrodynamic manoeuvring derivatives in way of the horizontal plane. Accordingly, hydrodynamic coefficients were evaluated by normalising ship design parameters against the characteristics of those available in the software database. It is noted that MCOL can use all hydrodynamic properties from any hydrodynamic simulations (3D, panel method, etc.). However, only restoring forces were calculated from 3D hull geometry.   For the grounding benchmark study, Fluid-Structure Interaction (FSI) was implemented by MCOL [37,38] for the pre-defined grounding impact speed defined in Table 4. The influence of resistance and damping forces was ignored. This is because they make the numerical idealisation uneconomic and influence minimally the results [22]. To calculate hydrodynamic properties, AALTO used the approach of [30]. BV/ICAM and MSRC utilised Hydrostar [33].  Table 8 summarise the applied extents and associated body boundary conditions. The participants modelled bulkheads for all scenarios. The majority of structural members of the target ships utilised in this study are made of mild steel. Limited amount of   structural elements above the main deck have been made of aluminium or anisotropic material (e.g., composite). On this basis it was assumed that mild steel is the main material for all structural members. Table 9 presents material properties and frictional coefficient. AALTO and MSRC used Cowper-Symonds equation. This formula is originally derived for the upper yield stress and accounts for the influence of dynamic effects on yield stress [39] (C = 40.4 and q = 5 for mild steel). All participants applied 0.1 of plastic strain as fracture criteria for all scenarios. BV/ICAM simulations did not take under account the friction. The hydrodynamic properties outlined in section 4 were applied as boundary conditions in the simulations of all participants. Both AALTO and MSRC used the LS-DYNA explicit code [40] for structural dynamic analysis. AALTO and MSRC modelled the plated structure by shell elements according to Belytschko-Tsay [40]. In the AALTO model stiffeners were modelled by 1D beam elements encompassing the Hughes-Liu cross section integration model [40]. On the other hand, for both striking and struck ships, MSRC modelled all stiffeners greater than 200 mm using shell elements. Following a convergence study, AALTO used 150 mm of element size in way of the expected damage area; see Kim et al. [22]. For the same analysis [3,23,41], MSRC selected elements size of 200 mm and 175 mm for striking and struck ships respectively. Their discretisation choice is based on engineering experience and reflects a quarter of frame spacing of the ship (e.g. 800 mm for the striking ship and 700 mm for the struck ship).

Shipship collision simulations
In the BV/ICAM method, both striking and stuck ships were modelled by very large-sized structural element units (plate, stiffener, girder, frame, etc.) and a limited number of nodal points (the so-called super-elements) in SHARP program [42,43]. Numerical calculations used experimentally validated closed-form analytical formulas of the resistance of each super-element for oblique collisions [44,45]. Then, by combining the individual resistance, it was possible to demonstrate the ability of Ships A, B to withstand a collision event. Examples of using SHARP super-element code in conjunction with MCOL solver can be found in papers by Paboeuf et al. [46,47]. In the super element method used by BV/ICAM, the fracture strain is not affected by the element size. This is because super elements for different structural members (plate, stiffener, girder, etc.) were developed based on experimental and analytical methods and accordingly each super element is linked with an empirical formula describing strain variation.
In a similar fashion MARIN modelled the complex ship structure as an aggregate of multiple constructional sub-elements the Fig. 7. Extent of analysis and structural layout in the grounding study.

Table 10
Dimensions of the analysis extents.  responses of which were previously defined analytically [34]. These sub-elements utilise closed form expressions derived by Refs. [44,48,49] and were used to calculate collision resistance.

Ship grounding simulations
In grounding simulations with forward speed, the high complexity of very sharp bow shapes and bow structures may link with geometric uncertainties. This is the reason why the simplified bottom structure after the bow of the target ship shown in Fig. 7 was used. Even if the structure for FE analysis is simplified, the ship's information of Ship A is used to calculate hydrodynamic properties. Table 10 summarises the extents of analysis including structural layout and details. All approaches presented accounted for the hydrodynamic properties outlined in section 4. Boundary conditions were implemented by MCOL. There were no body fixed boundary conditions. The assumptions for material properties and frictional coefficient used in grounding simulations are given in Table 11. AALTO and ICAM applied a bi-linear (elastic-or rigid-perfectly plastic) material curves. MSRC employed a true stress-strain curve based on experimental tests of Paik [50] suitably modified for the material properties shown in Table 11.
AALTO used element size-and plate thickness-dependent fracture strain for the definition of failure criteria (0.103 for inner bottom, 0.15 for longitudinal girder, and 0.126 for others) as per Eq. (2) [51].
where ε f is the fracture stain, t is the plate thickness, and l e is the medium (diagonal) element size, respectively [51]. Following a mesh convergence study element size of 80 mm was considered satisfactory. In ICAM simulations using solver FLAGS, only the bottom floors were associated to a failure criterion with strain threshold value equal to 0.2 which is validated with numerical and experimental methods. MSRC applied Eq. (3) to define the failure strain. This accounts for the influence of strain-rate effect on fracture strain. To calculate the strain rate in Cowper-Symonds equation [39], MSRC accounted for maximum strain rate during accidents. This is a function of the initial impact velocity defined by Eq. (3) [41]. The fracture strain values of 0.052 and 0.049 were applied for 5 knots and 10 knots, respectively. . 8. (a). Results of mesh convergence study by AALTO Fig. 8(b). Results of mesh convergence study by MSRC.
S.J. Kim et al.
where ε fd and ε fs are fracture strains under dynamic and static loadings, ξ is the ratio of the total energies to rupture for dynamic and static uniaxial loadings, ε is the strain rate, and V 0 is the initial impact velocity in m/s. respectively [39,41]. AALTO and MSRC adopted the LS-DYNA explicit code [40] with the same modelling techniques used in the collision analysis. Fig. 8 shows the results of their mesh convergence study for grounding scenario 1. Both participants used body fixed boundary conditions (side shells are fixed in all directions except for the surge direction) instead of MCOL to conclude on the most appropriate element size. AALTO selected 80 mm of element size, and MSRC used 250 mm of element size. In the explicit analysis of LS-DYNA, element erosion (fracture) that depends on the element length and plate thickness may lead to lack of convergence in the numerical solution when the same element fracture criteria are used. Accordingly, the element size-and plate thickness-dependent fracture strains calculated by Eq.
In ICAM's method, grounding dynamics were modelled by super-elements accounting for closed-form analytical formulations of the resistance of each unit. This is similar to the modelling approach used in collision simulations. Then, by combining the individual resistances, the dimensions of the breach were evaluated. The principles of FLAGS super-element grounding solver are summarised in  the simplified chart presented in Fig. 9. It is noted that the super-elements involved in ship bottom grounding simulations were idealised as solid plates between support members. Thus convergence study was not required. Depending on rock shapes (e.g., conical, paraboloid, etc.), FLAGS adopts different analytical expressions or closed-form solutions developed by Simonsen [52] and Sun et al. [53] for conical shaped rock, and Pineau and Le Sourne [54] and Pineau et al. [55] for paraboloid shaped rock.

Collisions
AALTO_2 demonstrate additional collision simulations which consider only the surge velocity of striking ship (see Table 7). Results show the effect of initial impact velocity as compared to 'AALTO_1'. Fig. 10 demonstrates a comparison of collision-induced penetration and internal energy. From an overall perspective, the penetration in oblique collision is smaller than perpendicular collision. This could be attributed to the effects of friction and the FEA modelling idealisation. The internal energy accounts for both absorbed energy from deformed and fractured structures and sliding energy due to friction and contact. For each of the super-elements implemented in SHARP solver, virtual displacement was implemented by experimentally validated closed-form expressions accounting for structural resistance (for example the plies and the plastic hinges formed in a concertina deformation mechanism of a deck or a bulkhead) [42][43][44][45]. Thus, the model embodied both plastic deformation and friction mechanisms. This is the reason why SHARP solver does not differentiate which part of impactor kinetic energy is absorbed by mechanical deformation and which part is dissipated by friction with the impactor or by self-contact.
Figs. 11-14 summarise results from the collision benchmark. Penetration values and energy differences are similar between BV/ ICAM and AALTO. Yet, in the AALTO method, it seems that the exact balance between mechanical and frictional energies as given by LS-DYNA output is quite sensitive with respect to contact modelling, the plastic behavior of structure and friction coefficient. As the total energy given by ICAM should represent both mechanical and frictional dissipation comparisons should account for both internal and frictional energies. The deviations in the results under scenario 3 relate with frictional mechanism. Also, the effect of body condition of the striking ship (deformable or rigid) is minor as the bow structure of striking ship is very sharp. Therefore the energy dissipation by structural response of deformable striking ship is smaller than the absorbed energy on the struck ship [22].
MSRC results show good agreement for perpendicular collision cases at low speed (scenarios 1 and 4) as compared to other participants. On the other hand, the MSRC results demonstrate deeper penetration at high speed collision (scenario 2) and higher internal energy for oblique collision (scenario 3). Structural crashworthiness results by MARIN and ICAM are in good agreement for scenarios  Figs. 11 and 12). This is because both software programmes use similar super-element techniques. However, for collision scenario 3 (oblique collision case), internal energy and penetration values predicted by MARIN are lower than those of ICAM. For most scenarios the two super-methods give close reaction forces. An exception to this is the oblique collision scenario during which MARIN's approach doesn't account for the influence of horizontal in-plane effects. Additional differences in the results displayed for scenario 4 could be attributed to the incomplete definition of transverse bulkhead super-elements implemented in MARIN code MarcolXMF. These observations are also justified by comparing MARIN against AALTO LS-DYNA simulations.

1,2 (see
Comparison of AALTOs results shows that the influence of 6 DoF initial velocity is minor because the struck ship is stationary in this benchmark study. Table 12 compares the timing of the first rupture of struck ship's side shell obtained by AALTO_2, BV/ICAM and MSRC simulations. All these model use the same hydrodynamic properties. However, the results show much different values of timing, penetration, internal energy and reaction force. This could be attributed to different FE idealisations. The super-element methods by BV/ICAM and MARIN, give similar results for scenarios 1 and 2. However, for scenarios 3 and 4, BV/ICAM and MARIN models display different mechanism for rupture initiation. This could be attributed to differences in ruptured structural members in oblique collisions (see side shell and the upper deck rupture sequences in BV/ICAM and MARIN approaches in scenario 4).  Table 13 compares the extents of damage by AALTO and BV/ICAM methods which represent nonlinear FEM and super-element method. It shows that the super element method gives similar damage extents in both penetration and horizontal damage. It is noted that when the collision speed is slower, the 6 DoF initial impact velocity leads to smaller damage extends. The influence of initial velocity in 6 DoF appears to be minor.  In general, results by ICAM in Figs. 17 and 19 appear to be in very good agreement with LS-DYNA simulations by AALTO for the case of the paraboloid rock. It appears that ICAM's method describes well the real material behavior (true stress-strain curve) for the grounding with paraboloid rock (see comparisons of ICAM_2 and MSRC_2), even if the rigid-perfectly plastic material is applied in the solver. However, there is a small difference between results by AALTO and ICAM_1 when the ship is grounded with the conical rock. In conclusion it appears that different hydrodynamic properties (e.g., the added mass in MCOL by AALTO is 10 times higher than the one in ICAM's MAOL) may lead to different structural behavior and the effects of FSI introduce additional modelling uncertainties. Table 14 shows a comparison of maximum damage extents between AALTO and ICAM simulations. It shows that the uncertainties from the calculation of hydrodynamic properties because ICAM_2 with MCOL input by AALTO gives closer results with AALTO than ICAM_1 which uses hydrodynamic properties by ICAM.

Table 12
The timing of the first rupture of struck ship's side shell after the first contact (T: time, P: penetration, IE: internal energy, RF: reaction force). Scenario    The approach by AALTO yields a smaller grounding length and higher internal energy than MSRC_1 due to the smaller fracture strain used in this model (see Figs. [16][17][18][19][20]. However, when MSRC_2 used MCOL by AALTO and fracture criteria by Eq. (3) (same approach with AALTO), grounding length reduced and the internal energy rapidly increased in short time. It may therefore be concluded that when similar fracture criteria are used, results by all approaches are in a good agreement for all scenarios. This demonstrates that the most influential factor on the difference of results by MSRC and others is the fracture strain. Thus longer grounding length and smaller internal energy is evident in MSRC_1 in comparison to AALTO simulations. On the other hand, when MSRC uses same fracture criteria with AALTO, the results are closer.
Numerical comparisons demonstrated that only 4 min are necessary to simulate grounding scenario 5 (ship running aground at 10 knots) using FLAGS/MCOL. On the other hand, 4 days are required when using parallel computing (12 of CPUs and 20 GB of RAM) for the FEA LS-DYNA/MCOL simulations.

Conclusions
This paper presented shipship collision and grounding benchmarks for structural response analysis of ships. The modelling methods accounted for hydrodynamic actions associated with operational and environmental conditions. The benchmark highlighted the importance of suitably justifying modelling assumptions to reduce uncertainties in modelling and simulation. Conclusions can be summarised as follows:  • In general, all approaches give quite similar maximum response in both collisions and groundings except for oblique collision and groundings accounting for sharp rock. However, when all participants use a same MCOL input and similar fracture criteria (BV/ ICAM_2 and MSRC_2), results of grounding scenarios 2 and 4 are also in a good agreement. This implies accurate modelling of fluid structure interactions is essential. • The extend of FEA model used for the analysis (i.e. partial or full ship FEM) may lead to big differences in terms of collision penetration. This is because of the combined effects of ship dynamics; i.e. the struck ship experiences both large sway and yaw motions while the striking ship is subject to surge and yaw motions. • Fracture strain critically affects structural responses. In general, the approach by MSRC_1 gives longer grounding length, but smaller internal energy than AALTO. This may be due to the smaller fracture strain and hydrodynamic properties (MCOL) in MSRC_1 assumption. However, when MSRC_2 uses Eq. (2) for fracture criteria as AALTO, for which the value is increased by 34% compared to MSRC_1, the grounding length becomes smaller. In this case the internal energy also rapidly increases in shorter time because the structure absorbs more energy by plastic deformation before rupture. Notwithstanding this further validation by largescale experiments for both collision or grounding cases would be beneficial. • In grounding, the different trend of results presented by MSRC for scenarios 1 and 3 (corresponding to conical shaped rocks) and scenarios 2 and 4 (corresponding to paraboloid shaped rocks) demonstrates the importance of using appropriate mesh sizing. This is also supported by the observation that a mesh convergence is essential to capture the crushing mechanisms during grounding. • The super-element solvers SHARP and FLAGS were shown to be accurate enough and efficient to solve collision and grounding mechanics of ships. They give quite similar damage extents in both accidents with nonlinear FE analysis. Accordingly they could be used for rule making within the context of SOLAS probabilistic ship damage stability requirements, where a large number of scenarios have to be investigated [29,56].

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.