A two-way coupled FSI model for the rapid evaluation of accidental loads following ship hard grounding

The assessment of structural crashworthiness following ship grounding events should consider the influence of multi-physics on dynamic response. To date, this has been achieved by computationally expensive approaches combining explicit 3D FEA with hydrodynamic solvers (e.g. Conti et al. (2021); Kim et al. (2021)). To offer an alternative, this paper presents a rapid six degree of freedom (6-DoF) two-way coupled fluid– structure interaction (FSI) model that could be used for the rapid evaluation of ship dynamic response and structural crashworthiness during hard grounding events. The time domain solution is based on a simplified contact analysis model that combines the external dynamics idealization of Taimuri et al. (2020) with an extended version of the internal mechanics analytical formulation of Simonsen (1997a) and Sun et al. (2017). The plate cutting angle is calibrated based on the rock geometry and hull indentation. A ray-tracing algorithm that utilizes panels and the tip of a conical rock is implemented to idealize hull penetration. The model is numerically validated by comparing results against LS-DYNA simulations for a box-shaped barge of double bottom configuration and a passenger vessel. Parametric studies demonstrate that longitudinal and vertical forces produced by different models compare well. However, the lateral forces for the study case of a passenger ship progressing along a straight path are shown to be underestimated. This is because the proposed method (i) does not allow for the influence of lateral loading on longitudinal hull members and (ii) overestimates the forces in regions of extreme curvature (e.g. fore shoulder, bow and boundaries of the hull).


Introduction
Ship collisions and groundings are of great significance as ship accidents persist regardless of regulatory progress. Allianz Global (2020) reports that ship groundings have been the second biggest cause of ship losses during the past decade. Factors contributing to ship grounding accidents relate to machinery failures, human factors, traffic complexity and environmental conditions. Consequences may be foundering leading to sinking or submerging of marine vessels, environmental pollution, traffic flow interruption and loss of human life.
Grounding events arise when ships crash into the seabed or run onto protruding rocks. Hence, ship grounding dynamics are complex and interlink with maneuvering forces, environmental loads emerging from short waves, wind, ocean currents and an increase of ship resistance. The dynamics of ship grounding can be studied by model tests or fully coupled fluidstructure interaction (FSI) methods. The former are very expensive, time consuming and may lead to results that are  TYc,global F TZc,global Total contact force in a global coordinate system (N) F V Vertical resistance and reaction force of structure (N) Total number of block in x, y, and z direction for storing panel ids Transverse and longitudinal initial metacentric height. g Acceleration of gravity (9.807 m/s 2 ); horizontal to plastic force ratio H def Height of the transverse floor deformation (m) I X , I Y .I Z Roll, pitch, and yaw moment of inertia (kg m 2 ) K Restoring force matrix K G Vessels' CoG measured from keel (m) Kṗ Variation of roll moment due to roll acceleration (kg.m 2 ) Kṙ Variation of roll moment due to yaw acceleration (kg·m 2 ) Kv Variation of roll moment due to sway acceleration (kg·m) k vertical to horizontal force ratio loc floor Location of the transverse floor from base of the rock (m) L pp Length between perpendiculars (m) M

Inertia matrix M A
Added inertia matrix M 0 Plastic bending moment of plate per unit length (Nm) Mq Variation in pitch moment due to pitch acceleration (kg·m 2 ) m Mass of the ship (kg) mod computes the remainder of the division N 0 Plastic membrane force of plate per unit length (N) Nv Variation in yaw moment due to sway acceleration (kg·m) Nṙ Variation in yaw moment due to yaw acceleration (kg·m 2 ) Location of center of flotation measured from CoG (m), +ve bow x G Longitudinal distance of the center of gravity measured from midship (m), +ve bow (m) x int , y int , z int

Location of rock tip from ship origin (m) Y Hull
External sway damping forces (N) Yṗ Variation in lateral force due to roll acceleration (kg·m) Yṙ Variation in lateral force due to yaw acceleration (kg·m) Yv Variation in lateral force due to sway acceleration (kg) Effect of sway and yaw linear velocities on lateral force Y ′ vrr , Y ′ vvr Coupled effect of sway and yaw velocity on sway force and yaw moment Y ′ vvv , Y ′ rrr Changes in lateral force (Y ′ ) and yaw moment (N ′ ) due to higher-order pure sway velocity and yaw velocity.

Zq
Variation in heave force due to pitch acceleration (kg·m) Zẇ Variation in heave force due to heave acceleration (kg) z G Vertical distance of the center of gravity measured from midship and still waterline +ve towards  Kitamura and Kuriowa, 1996;Lemmen et al., 1996;Rodd, 1996a,b;Rodd and Sikora, 1995). In theory emerging multiphysics methods can provide a useful alternative (e.g. Rudan et al., 2019;Song et al., 2016;Kim et al., 2021). However, with the exception of the Super Element method they are computationally uneconomic (Kim et al., 2022b). The latter is important, especially considering that the mitigation of grounding risks requires accurate and fast estimation of hull damage extents and the development of Risk Control Options for use in ship design for safety (Conti et al., 2021;Zhang et al., 2021).
To offer a scientifically sound yet practical alternative this paper introduces a rapid time-domain two-way coupled FSI model. The method combines 6-DoF ship maneuvering motions with a model accounting for structural deformations and the influence of hydrodynamic properties in way of grounding contact. The assessment of structural deformation is based on the energy dissipation approaches of Simonsen (1997a) and Sun et al. (2017). The multiphysics crashworthiness model presented advances the state of the art by allowing for: • The influence of plate split angle variations at different penetration depths. Consequently, the splitting angle of the plate is defined as a function of rock penetration and base-radii. Likewise, the extent of plate deformations depends on out-of-plane ship motions.
• The variation of the inner and outer bottom plate split angles that were considered equivalent in past models.
• The evaluation of dynamic rock-ship interactions, through a rapid coupling algorithm that combines internal and external mechanics at the interface of the rock tip with the ship's bottom.
Results validated for the case of a passenger vessel subject to hard grounding with a conical rock suggest that the method is computationally robust, economic while it accounts for environmental loads (wind, shortwaves, and ocean currents). Thus it could be used for low-fidelity ship design and optimization, to introduce new damage extend datasets, rules and regulations for ship damage stability assessment and for the development of intelligent decision support systems for use in emergency response.
The paper is organized in five sections as follows. State of the art research in ship grounding is reviewed in Section 2. The fundamentals of the novel FSI method and associated assumptions are explained in Section 3. The validation of the rapid FSI model against LS-DYNA simulations and discussion of the results are included in Section 4. Conclusions are drawn in Section 5.

Literature survey
A review of several applications implementing FSI methods for the evaluation of loads and responses is given by Hirdaris et al. (2014) and Hirdaris and Mikkola (2021). The evaluation of sea loads associated with accidental events such as ship collisions (e.g. Choung et al., 2020;Lee et al., 2017;Rudan et al., 2019;Zhang and Suzuki, 2006), groundings (e.g. Kim et al., 2022bKim et al., , 2021Kim et al., , 2020Lee et al., 2017) and ship-ice interactions (Lu and Zou, 2020;Luo et al., 2020;Song et al., 2016;Wang and Derradji-Aouat, 2010) account for the influence of contact mechanics by novel FSI models. The use of such multiphysics methods for the assessment of ship loads in intact or damaged conditions requires considerable computational time because of two-way coupling and the fidelity of FEA models Lakshmynarayanana and Hirdaris, 2020). This is the reason why the derivation of simplified modeling approaches (e.g. empirical methods, super element methods, and contact mechanics models) remains relevant within the context of developing advanced regulations for ship design assessment and the risk management of shipping operations (e.g. Idrissova et al., 2019;Descantes et al., 2019;Conti et al., 2021;Zhang et al., 2021).
Modeling grounding dynamics consists of two processes namely external dynamics and internal mechanics. A summary of published work on the subject is given in Table 1. Accounting for the influence of hydrodynamics in way of contact (e.g., added mass and damping effects due to ship motions, evasive actions following a maneuvering path) relates to external dynamics. Internal mechanics mostly encompass structural responses and deformation because of energy dissipation post-mortem a grounding event. The influence of evasiveness has been traditionally ignored in these models. The latter may be considered disadvantageous especially considering that as a ship maneuvers in shallow waters a wide range of seabed topologies namely rocks, shoals, and reefs may be evident (Alsos and Amdahl, 2007). In literature, rocks encountered in hard grounding scenarios are idealized as pinnacles or truncated cones that will rupture the hull bottom. Shoals have large and soft contact area (Ehlers, 2011) where typically a hull girder may collapse and the bottom structure may experience crushing and denting. Reefs can be idealized as conical rocks (Simonsen, 1997b). A ship's striking on a rock, or a shoal is termed power grounding (bottom raking and bottom sliding). In some occasions, stranding may be followed by grounding (Rizzuto et al., 2018;Souliotis and Samuelides, 2013).
As part of modeling internal mechanics a ship's structural deformation can be predicted by empirical methods (Minorsky, 1958), FEA (Kim et al., 2020;Liu et al., 2021), the super element method (Buldgen et al., 2012;Kim et al., 2022b;Le Sourne et al., 2012), experimental approaches (e.g. Calle et al., 2020;Rodd, 1996a,b;Rodd and Sikora, 1995) or simplified analytical models (Zhang et al., 2019). External dynamics have been modeled by the so-called decoupled methods that traditionally tend to ignore the influence of fluid dynamics and predict the energy due to impact by analytical expressions or FEA (e.g. Liu et al., 2017;Yu et al., 2019;Ehlers, 2011;Pill and Tabri, 2011). Decoupled FEA simulations have been used for ship structural design optimization (e.g. Song and Hu, 2017;Zeng et al., 2016), to idealize plastic deformation (Abubakar and Dow, 2013) and to introduce empirical formulae emerging from force penetration curves (Heinvee and Tabri, 2015). Recently, Kim et al. (2020Kim et al. ( , 2021 presented large scale numerical crashworthiness analysis methods that account for the  (Rodd, 1996a,b), Pan: Panel model of full ship without structure details; HG: hard grounding; SG: Soft Grounding.; Dry: without external mechanics. effects of ship grounding dynamics. This work demonstrates that the influence of hydrodynamic restoring forces may also be significant and if coupling effects are ignored the ship may ascend erratically after her initial contact with the seabed. The computational expense of the above mentioned large scale multiphysics methods, implies the need to develop simplified approaches for use in rapid design development and assessment. Lately such reduced order integration methods have been presented by Kim et al. (2022a,b). The former makes use of the so-called Super-Element method of Le  according to which a structure is divided into several macro elements (e.g., plates, intersections, transverse and longitudinal members) and forces associated with the structural resistance of each of those is described by closed-form analytical expressions. As elements deform individually, the aggregate of their deformation is used to express the total structural deformation of the ship in accidental conditions. The paper by Kim et al. (2022a) demonstrates an alternative approach that combines simplified spring elements for the idealization of restoring forces and an FEA model for the idealization of grounding dynamics.
A critical review of the above leads to the conclusion that relatively few studies demonstrate the influence of direct coupling of mechanics while accounting for hydrodynamic actions. An exception to the rule are the recent models presented by Kim et al. (2022a,b) and Lee et al. (2017Lee et al. ( , 2013. The papers by Lee et al. (2017Lee et al. ( , 2013 and Kim et al. (2021) demonstrate that strongly coupled ship grounding dynamics can be idealized by advanced discretization solvers such as Arbitrary Lagrangian-Eulerian (ALE) and nonlinear FEM possibly with spring element corrections to account for the influence of restoring forces. Yet computational economy is an issue with such type of large-scale models. On the other hand, the work by Kim et al. (2022b) confirms the emergence of the Super Element method as a possible alternative. However, even this approach should be validated further and maybe inevitably inaccurate in terms of modeling complex geometries, rock dynamics and ship dynamics.

FSI method
In this work, the 6-DoF dynamics of a rigid ship are coupled with an analytical contact mechanics model. Coupling is implemented by a contact detection algorithm. An enhanced technique is proposed to determine the instantaneous change of plate splitting angle. The model has been validated by LS-DYNA FEA simulations encompassing a hydrodynamic model that accounts for the influence of hydrodynamic inertia, damping, memory effects and restoring forces and moments via the 'Mitsubishi Collision -MCOL' (Ferry et al., 2002;Le Sourne, 2007) algorithm. Fig. 1 illustrates the process of solving a simplified two-way coupled partitioned approach that transfers information from external to internal mechanics at each time step. Ship maneuvering data, structural information and rock parameters are accounted for by the method. At the initialization stage, the shape of the hull is used to form 3D sub-blocks and store the panels. The ship maneuvers with controlling devices. Thus, interaction with the ground is detected based on information available for the rock tip location and ship position. Ship motions are encompassed in the coupling algorithm. In way of an FSI interface the rock penetration (displacement) inside the hull can be determined and the structural deformation is used to compute the contact forces. The structural solver makes use of the momentum equation to estimate the contact forces and moments at the FSI interface. Ship motions are updated at each time step.
External dynamics are solved using the approach of Taimuri et al. (2020). Accordingly, the ship can be modeled as a rigid structure maneuvering in deep or shallow waters, in calm or short-wave conditions and with or without wind and ocean currents. Shallow waters have a larger influence when the water depth to ship draft ratio is between 1.5 and 1.2. A ratio of the order of 3 (or higher) corresponds to deep water conditions (Vantorre et al., 2017). According to HELCOM -Helsinki Commission (2021) grounding statistics, sea depths range from 10 to 500 m in areas where grounding occurred. The grounding accidents of Exxon Valdez and Costa Concordia indicate that both deep and shallow sea conditions should be considered when assessing grounding risk. Notwithstanding this, the present study considers only the case of a deep-water scenario.
FSI is governed by the rock's tip and the triangular hull panels. To idealize the physics of plate tearing the hull panels were stored in identical sub-blocks specified by 3D Cartesian coordinates that bound the complete hull shape. The required sub-block can be identified directly from the coordinates of the rock point on a ship-fixed coordinate system. For those cases that the rock is located within a sub-block, a ray tracing algorithm can be used to detect the intersection of a vector heading towards the ship base (Moller and Trumbore, 1998). This modeling strategy helps to accurately estimate the distance between the hull panel and the rock tip, which may then be used to calculate contact forces. Once an intersection is detected, the structure will deform. Thus, internal mechanics make use of the so called ''upper-bound theorem of plasticity'' (Simonsen, 1997a;Sun et al., 2017). Output results from this numerical scheme are the ship trajectory, grounding length, absorbed energy by the ship and forces leading to deformation. Different rock geometries can be formulated based on the gap openings and flap angles of different shapes (Simonsen, 1997c,b). However, the method presented in this paper is only applicable to conical rock shapes with rounded tips.

FSI modeling assumptions
The internal mechanics model (Simonsen, 1997a;Sun et al., 2017) introduced in this paper is based on the following assumptions: • The structural material is rigid and plastic. Accordingly, elastic deformation is neglected and von Mises yield criteria are followed; • The rock is idealized as a cone with a rounded tip radius and a semi-apex angle; • In the first instance, the deformation mechanism of each of the structural components is considered separately.
Following this, the resistance of all structural members within the deformation zone is summed up; • Before fracture, the plate follows a prescribed deformation mode according to which the plate wraps around the rounded tip of the cone; • After fracture, the bottom structure of the plate conforms to the rock and the generatrix (line segment between cone base and apex of the cone) is used to define the contact force on the plate fractured; • The plate cutting phenomena follow stable or clean curling cut. Concertina tearing and braided cut of the plate are not considered; • Strain effects in the longitudinal direction of the plate, girder and stiffener components are not accounted for; • The plate plastic resistance and transverse floors are modeled according to Simonsen (1997a), Sun et al. (2017) respectively. Accordingly, (i) longitudinal girders and bulkheads are unstiffened and only their membrane effects are considered; (ii) longitudinal stiffeners only account for bending deformation; • The deformation zone ahead of the rock is defined by penetration into the hull and a cone semi-apex angle.
According to Taimuri et al. (2020) assumptions related to external mechanics are summarized as follows: • In-plane hull hydrodynamic damping coefficients are based on semi-empirical expression that originates from the regression of captive model tests (CMT). Out-of-plane radiation damping is estimated using linear theory (natural period of the ship).
• Added inertia can be approximated either by a semi-empirical expression (Brix, 1993;Clarke et al., 1982) or strip theory (Frank, 1967;Journée, 1992); • Shallow water effects are accounted for by an increase in resistance. However, the influence of wave-making resistance, squat and bank effects are not considered; • Added inertia and damping in shallow waters are modeled according to Vantorre (2001) Heave and pitch motions are not considered; • Short waves account for mean 2nd order wave drift forces (Faltinsen, 1990;Faltinsen et al., 1980;Sakamoto and Baba, 1986). The effect of limiting water depth on the 2nd order mean drift force is defined according to Longuet-Higgins (1977). LS-DYNA simulations assumed that all structural components are made of mild steel (S235). Plates were modeled by 6 -DoF shell elements with Belytschko-Tsay formulation (LSTC, 2021) and the Hughes-Liu degenerative continuum model was used to model stiffeners, girders and frames (LSTC, 2021). Based on the convergence study of Kim et al. (2022bKim et al. ( , 2021, element sizes of the order of 150 mm and 100 mm were considered satisfactory in terms of idealizing the bottom plate tearing for the passenger ship and the box barge models respectively. The Cowper-Symonds equation was applied to model strain effects (Jones, 2011). As part of this process, uncertainties of relevance to the estimation of the dynamic fracture strain factor have been thoroughly considered. According to Paik et al. (2017), the dynamic fracture strain rate for mild steel is around 0.07 for 14.5 mm plate thickness and 191 mm of element size. Thus, for the topology considered in this paper the fracture strain factor would be around 0.087-0.093. On the other hand, Lehmann et al. (2001) empirical formulae suggests a fracture strain of 0.11-0.12. In conclusion, an average fracture strain of the order of 0.1 was considered adequate .
The external mechanics model used in LS-DYNA-MCOL coupling simulates the case of a 6-DoF rigid ship heading straight with an initial forward speed in deep and calm waters. The added mass has been considered to remain constant at infinite frequency. Hydrodynamic radiation damping with memory effects were accounted for. The simplified 6-DoF maneuvering model used the same (infinite frequency) added inertia coefficient and ship heading with forward speed. However, it did not account for the influence of controlling devices and resistance effects. In-plane hull hydrodynamic forces were modeled according to Taimuri et al. (2020) and out-of-plane linear damping was evaluated as described in Appendix A.

External mechanics
The coordinate systems used for (i) the idealization of ship dynamics and (ii) the direction of structural deformation forces due to the penetration of the rock into the hull are shown in Fig. 2. Ship motions were defined using the Earth fixed inertial coordinate system ⃗ R = [X, Y , Z ] associated with the navigational position of the vessel. The body-fixed coordinate system (r 0 = {x, y, z}) was attached to the intersection of the plane passing through the waterplane area amidships the symmetric hull . The rock was specified as a point node, in the inertia coordinate The governing equation of ship dynamics (Fossen, 2011;Matusiak, 2021;Taimuri et al., 2020), structural deformation (Simonsen, 1997a;Sun et al., 2017) and contact detection is: where,Ẍ = {u,v,ẇ,ṗ,q,ṙ} T ,Ẋ = {u, v, w, p, q, r} T and X = {x, y, z, φ, θ , ψ} T are acceleration, velocity, and position vectors (translational and angular) in the body-fixed coordinate system; M and M A are the inertia and added inertia matrix respectively; C is rigid-body Coriolis and centripetal forces and moment matrix showing the rotation of the body in the inertial frame as the ship's origin r 0 is at a distance ⃗ r G = [x G , y G , z G ] away from the center of gravity; B P is the radiation damping matrix, B H is hydrodynamic linear and higher-order coefficients considering viscous damping, K is the restoring forces and moment matrix; F prop , F Rud , and F sw are the contributions of propeller, rudder, and short waves respectively. The contact forces and moments are defined by the vector F C . The details of inertia, Coriolis, centripetal, restoring forces and the damping matrix are given in Appendix A.
Structural deformations were coupled with 6-DoF ship external dynamics model via a contact detection algorithm (see Section 3.4). The governing 2nd order non-linear differential Eq. (1) was numerically solved using the 4th order Runge-Kutta time-stepping scheme. The response of the ship was calculated in the body-fixed coordinate system. To transfer the response into the global/Earth-fixed coordinate system the transformation matrix T (γ) shown in Eq. (2), was utilized: where γ represents Euler angles, T T (γ) is the transformation matrix of translational velocities from body coordinate system to inertial coordinates and T R (γ) is a transformation matrix of angular velocities into the inertial coordinate system. In shorter notation Eq.

Internal mechanics
The contact force F C (see Eq. (1)) results from impact loads. In the models of Simonsen (1997a) and Sun et al. (2017) the energy dissipation is estimated by integrating the stress and strain field over the volume of the structure to form an analytical equation. The structural resistance force presented in Eq. (4) reflects the energy balance in way of grounding contact; i.e. it expresses the relationship that the rate at which energy dissipates internally by plastic deformation, friction, and by fracture must be equal to the external work done.
where, F H is the resisting force of the structure in the direction of U; U is the relative velocity between the ship and a rock;Ė p ,Ė c ,Ė f express the rate of plastic, crack, and frictional energy dissipation respectively; F P is the plastic resistance which includes both plasticity and fracture; µ is the Coulomb coefficient of friction; The last term of Eq. (4), represents the sliding energy dissipation due to normal and frictional forces in way of the contact surface of the rock and the plate (p is the plate's normal pressure on the rock based on contact area S between rock and plate; U rel is the relative velocity between rock and plate according to Simonsen, 1997a). Hard grounding was assumed to take place over the rigid conical rock with a rounded tip radius R R and a semi-apex angle φ (see Fig. 3). The forces generated by the individual structural components were summed up to produce the resultant contact forces. The hull outer shell holds the structural components (floors, girders, bulkheads, and stiffeners) which are attached to the shell plate. Therefore, the global mode of deformation is governed by the outer hull plating, where the intersections are assumed to be intact during the deformation process. The generalized deformation modes and plate openings following grounding can be idealized as shown in Fig. 4. The membrane strain can be defined by the longitudinal and transversal gap openings u 0 and v 0 respectively (Simonsen, 1997a). A summary of the plate opening, extent of deformation, fracture penetration and structural component forces are presented in Appendices B and C. The horizontal and vertical reaction forces were evaluated by utilizing the plastic force, rock geometry and plate split angle. The effect of friction ''g'' and vertical to horizontal force ratio ''k'' were defined for both intact and fractured hull plating (Simonsen, 1997c). Eq. (5) shows the horizontal to plastic force ratio and vertical to horizontal force ratio estimation for the intact plate. Eq. (6) represents fractured plate.
where, µ is the coefficient of friction, ξ is the rolling of the plate and θ is the plate split angle.
The plate split angle represents the extent of plate deformation ahead of the conical rock. It helps demonstrate how far the length of the plate opening shall propagate before reaching a plastic zone. The plate split angle relies on structural component deformation energy and is critical in terms of suitably idealizing damage extents (see Eqs. (5)-(6)).

Plate split angle idealization
In Simonsen (1997a) the plate split angle θ (see Fig. 4) was defined as an unknown parameter that remains constant and minimizes the total resistance force of the structure. The plate split angle depends on the geometry of the indenter (rock). For instance, splitting of a plate can be assumed as the semi-angle of the wedge when using a wedge-shaped indenter. However, the splitting of a plate over a conical rock can vary depending on the depth of penetration and the actual shape. The following discussion highlights the importance of properly evaluating the plate split angle by demonstrating a comparison of plate split angle obtained from FEM solver against the theoretically minimized resistance force approach introduced by Simonsen (1997a). Fig. 5, illustrates the bare plate splitting angle approximated by the simplified (Simonsen, 1997a) model at different penetration depths with a cone semi-apex angle of 45 • . Fig. 6, shows the more detailed LS-DYNA FEA simulations with 2 m depth and cone semi-apex angle of 45 • . In the examined scenarios for LS-DYNA simulations the rock was assumed to be rigidly fixed in all directions and surge translations were applied at the side shells of the box (or plate). Plate specifications are presented in Appendix E. The bare plate split angle varied between 34 • to 48 • , depending on the time and position of the rock. The corresponding plate split angle at 2 m of rock penetration was around 6 • (see Fig. 5 encircled region). This value deviates quite significantly from FEM simulations and helps underline the importance of accurate split angle approximations.
To examine the effect of structural details on the variation of the split angle, a simple bare-plate and a complete boxship structure model comprising of double bottom, center girder, stiffeners, plate, and floors were simulated in LS-DYNA (see Fig. 7a, b) for the scenario presented in Appendix E. The original data was extracted at different positions along the length of the structure, where the velocity of the structure drops lengthwise. Simulations have been categorized by the base radius of the cone, height of the cone (penetration) and cone semi-apex angle. Fig. 7a, b show that the variation of a bare plate split angle is less sensitive to the speed drop of the plate in comparison to a fully structured ship. In a bare plate, the values are close to each other irrespective of cone geometry and penetration depth. Notwithstanding this, slightly more disperse data for a few categories were observed in full hull forms (e.g., δ out = 1 m corresponds to φ = 45 • , while δ out = 2 m corresponds to φ = 27 • ).
A summary of the results of LS-DYNA simulations is illustrated in Fig. 8. The minimum and maximum plate split for a bare plate are 26 • and 56 • respectively. For a complete box-structured ship these values become 32 • and 63 • (see Fig. 8). The mean variation of the plate rolling angle ranged from 0.56-0.72 times the plate split angle. These values can be used for estimating ξ . Theoretically the bound of the plate rolling ξ due to the relative velocity of the ship and contact surface  is in between 0 and θ (Simonsen, 1997a). For ease of use, in the present analysis plate rolling was assumed to be 0.5 times the plate splitting angle (ξ = 0.5θ ). These simulations underline the fact that, the approach proposed by Simonsen (1997a) which selects a plate split angle at minimum resistance force, leads to significant underestimation of the plate split angle as compared to FEA simulations.
The curve fitting of the plate split angle as a function of penetration into the hull and radius of the cone base is presented in Eq. (7) and Fig. 9. The curve fitting method used to determine the plate split angle has been formulated as a function of B de and δ out . As part of this process a complete ship model and average values of the plate split angle were where, a, b, c, and d are the coefficients of the polynomial equation, i subscripts correspond to the four different rock base radii. These four different curves were fit according to Eq. (8), as a function of the rock base radii B de (see Fig. 9 bottom).
Back substitution of the polynomial coefficients of Eq. (8) into Eq. (7) resulted in a split angle θ expressed as a function of B de and δ out . This new equation of plate split angle helped simulate the dynamic behavior of the ship structure as a function of conical rock geometry. The model is an improved version of the previously considered average value for the entire ship introduced by Simonsen (1997a) and, depending on penetration indentation, may lead to different horizontal and vertical forces (see Eq. (6)).

Contact detection algorithm
To ensure appropriate modeling of the penetration into the hull and rapid simulation of the grounding event, a fast contact detection algorithm has been developed. This algorithm produced the penetration depth from which the structural deformation was calculated at each time step. The implementation of the method is illustrated in Figs. 10, 11. The rigid hull model was idealized by triangular panels. The coordinates of each panel were defined from the ship's origin (see Fig. 2). The key steps of the algorithm are summarized as follows: Based on the ship's main dimensions a ship bounding cartesian grid that consists of identical 3D sub-blocks (Gx, Gy, and Gz along the length, beam, and depth of the ship respectively) is generated (cartesian grid initialization in (d) The rock tip is idealized as a point in an inertial coordinate system represented by a vector ⃗ R Gr (Fig. 2). The coordinates of this point are transformed into the body-fixed coordinate system, while the search of the rock within the sub-block is conducted in the ship fixed coordinate system; (e) The position vector from the ship origin to the rock tip is used to determine whether the rock is within the ship's bounding cartesian grid or not. If the rock tip is inside the cartesian grid, a 3D sub-block that contains the rock tip is identified Fig. 10. Otherwise, the solver returns to the equation of motion depicted in Fig. 11. (f) If the rock tip is located inside the 3D sub-block, the next step is to locate the intersection between panels and the tip of the rock. If there are no panels in a sub-block containing the rock tip, then neighboring blocks in the direction of the rock tip to the base are searched until an intersection with a sub-block grid is detected and stored. The depth of the penetration is calculated by a ray-tracing algorithm (Badouel, 1990;Moller and Trumbore, 1998).
(g) After the intersection is identified, the penetration into the hull δ out , is output from the contact coupling model.
This penetration is then transferred to the internal mechanics model to evaluate the forces leading to structural deformation.
The contact force on individual structural members in way of an intersection depends on the penetration and details of structural members at the location of the rock-ship interaction. The horizontal force F H is decomposed into the longitudinal and lateral components by using the relative velocity between rock and ship intersection panel. The velocity of the ship shell at the intersection point is calculated from the velocity vector field equation: where, ⃗ V int = [u int , v int , w int ] is the velocity of the panel at the intersection point; ⃗ V = [u, v, w] is the body velocity of the ship, ⃗ r int = [x int , y int , z int ] depicts the intersection position of the rock from ship origin and ⃗ ω = [p, q, r] is the angular velocity of the body. These velocities refer to the location of the rock and panel intersections, which are in a body-fixed coordinate system. The velocities are then transformed to the global coordinate system via the transformation matrix: The components of longitudinal, lateral, and vertical force components in the global coordinate system are defined as: F TZc,global = −F V In the above expression, the negative signs in the longitudinal and lateral planes express the relative velocity between the ship panel and rock. The total force is the sum of plastic plus fracture and frictional forces (see Eq. (4)). The breakdown of the total force into plastic resistance force (plasticity and fracture) and frictional force of the structural member is given is Eqs. (12)-(13). Plastic contact force, Frictional contact force, The origin of the plastic resistance force F P is demonstrated in Appendix B and the total horizontal and vertical forces are given in Eqs. (5) and (6) respectively. The total contact forces displayed by Eq. (11), in the body-fixed coordinates at the point of intersection, are expressed as: After the transformation of the forces to the body-fixed coordinate system, the contact forces become:

Input details for validation cases
The FSI model presented in Section 3 of this paper has been numerically validated via direct comparison with LS-DYNA-MCOL simulations. Results were compared against the constant plate split angle approach of Simonsen (1997a). Two study cases corresponding to (a) a simplified box-shaped ship and (b) the modern passenger vessel FLOODSTAND SHIP B (FSB) (Luhmann, 2009) were considered. The sub-cases of centerline and off-center (B/4) grounding were explicitly studied to better understand grounding dynamics under realistic scenarios. The plate split angle (present model and Simonsen approach), structural details, rock parameters, and ship initial conditions of box-shaped and passenger vessels are summarized in Table 2. Models of the box-like and passenger ships are shown in Fig. 12.
In the box-like ship model plate thickness and spacing of the floors were assumed to remain constant. Double bottom and uniform structural arrangements had only a centerline longitudinal girder. The passenger ship model consisted of more comprehensive and varying structural details namely different plate thicknesses, double bottom heights and location of the longitudinal girder (see Fig. 13). This variation was idealized in the simplified FSI model by dividing the ship into 6 longitudinal sections and by specifying structural details as illustrated in Fig. 13. The results presented in the following sections follow the coordinate setup adopted by LS DYNA as depicted in Fig. 12. The origin of the ship was assumed to lie on the center of gravity. The velocity transformation from the maneuvering model coordinates  to LS-DYNA coordinate system is presented in Appendix D. It is noted that the MCOL-LSDYNA coupling accounted only for heave, roll, and pitch restoring forces (see Eq. (5)). Therefore, in the validation study the effects of surge, sway and yaw motions in the restoring forces matrix (K) of Appendix A were ignored. The input details for MCOL simulations are given in Appendix F.

Numerical validation of a box-shaped ship centerline grounding
The time history of ship motions, sliding forces, and deformation energy of the centerline grounding case are summarized in Figs. 14-16. The top three plots in Fig. 14, represent the surge, heave, and pitch velocities of the ship. The bottom three figures depict their respective displacements. The longitudinal damage extents evaluated when using the present simplified FSI method were 4.8% shorter than those evaluated by LS-DYNA simulations. This resulted in marginally different damage lengths (36.80 m for FSI versus 38.65 m for LS-DYNA). On the other hand, surge velocity trends predicted by the simplified FSI and LS-DYNA simulation matched well. The slight shift of the slope in surge velocity after 2 s (see Fig. 14) reflects that as the ship is moving up, the penetration into the hull reduces, eventually leading to reduction of the longitudinal resistance force of the structure (see Fig. 15). Subsequently, the surge velocity slope changes when the ship starts to heave down after 4 s, and the longitudinal and vertical structural resistance increases. The values for heave and pitch motions of the simplified FSI model are higher in comparison to LS-DYNA simulation. This is because of the higher (i.e., more conservative) estimation of the forces by the simplified structural deformation model. At around 0.4 s the ship experiences the largest vertical force of 4.91 MN which leads to a higher heave response (see Fig. 14).
On average, the differences of the horizontal and vertical sliding forces predicted by the simplified FSI model and LS-DYNA are 0.6% and 20.0% respectively. At the outset of the contact, there is an abrupt jump of the forces predicted by the simplified model (see Fig. 15). In a real scenario, once the base of the rock is in contact with the hull, it progressively moves ahead to deform the structure. However, in the simplified approach introduced here, the contact is modeled based on the tip of the rock and the initial contribution of the rock base is ignored. This is explained in Appendix G. The rupture    Fig. 9. Eqs. (7) and (8) deg Simonsen approach plate split angle θ 11 10 deg of the transverse floor is seen as a sudden drop after each peak in the sliding force plot (see Fig. 15). The time to rupture increases as the ship velocity reduces and gets closer to zero. This behavior becomes more visible as an elongation in the force-time curve after 8 s (see Fig. 15). In the rapid FSI method proposed the deformation energy of the structure did not take into consideration the elastic behavior of the structure. Fig. 15, presents two internal energy components namely: (a) the sliding energy that comprises normal and tangential components of the contact force (frictional force), and (b) the total energy which is the sum of absorbing and sliding energy components 1 . The former is higher from the one of the simplified model and results in a 13.7% increase in results as compared to LS-DYNA. The total deformation energies predicted from both models show a rather marginal variation (less than 1%).
When the Simonsen concept of constant plate split angle is employed, the ship's longitudinal damage length increases significantly (Fig. 14). This is because the specified plate split angle lowers the total energy dissipation (Fig. 15). It reduces resistance to ship structure during grounding and prolongs the time for a ship to come to rest. The vertical sliding force obtained from Simonsen's approach (Fig. 15) 1 is slightly larger than the simplified model (Section 3.3.1). However, the overall energy decreases. This implies that the plastic resistance force is underestimated, and the frictional force is overestimated. Table 3, compares the damage extents predicted on the basis of the simplified model and Simonsen's approach against LS-DYNA simulation. For all cases the maximum damage width is located at the beginning of the contact when the vertical force is maximum. The minimum width is obtained when the ship has the largest heave. Moreover, the inner bottom in all cases remains intact. Fig. 16(a), demonstrates the damaged hull profile that resulted from LS-DYNA simulations. Fig. 16(b) and (c), illustrate respectively the damage width, length, penetration, and resultant forces from the present method and Simonsen's approach (1997a) respectively. Fig. 16(a) and (b) show good correlation of damage length and width between the proposed model and LS-DYNA simulations. However, 61% higher estimation of the damage length is evident when choosing the method of constant plate split angle (see Fig. 16(c) and Table 3).

Numerical validation of a box-shaped ship subject to off-center (B/4) grounding
The results of the off-center grounding (quarter breadth of the ship B/4) are shown in Figs. 17-19. The total damage length of the ship predicted by the simplified FSI model is 49.7 m. This overestimates the LS-DYNA predictions by 4.7%. The trends of the surge, heave and pitch motions are similar for both models. During off-center grounding alongside heave and pitch, the ship also experiences roll motion. For the box-like idealization presented here, the simplified FSI model predicts a maximum roll velocity of 3.5 deg/s and heeling angle of the order of 9 deg. The latter exceeds almost 1.25 times predictions by LS-DYNA. The average longitudinal, lateral and vertical sliding force variations predicted by the simplified model differ by 7.8%, 10.0%, and 11.2% respectively to those predicted by LS-DYNA. The increased sway motion in the simplified model is predominantly caused by the large yaw and somewhat influenced by the roll inclination. The maximum total energy and sliding energy differences lie within 0.5% and 10% to those predicted by LS-DYNA.
A graphical illustration of damage extents from LS-DYNA simulations against Simonsen and the new model are shown in Fig. 19. Table 4 summarizes the differences in damage extents. The simplified model shows that in the transverse direction the damage reaches the centerline longitudinal girder, as the sway displacement increases by 0.5 m after 7 s. This becomes visible from the extent of the damage width in way of 31 m length of the damage (see Fig. 19). An offcenter grounding demonstrates a similar tendency when a constant plate tearing angle (i.e., Simonsen's approach) is used. That is the substantial increase in damage extent may relate to a decrease in the total deformation energy (see Figs. 18,19(c)). Irrespective of the model used noticeably the inner bottom remains intact because the height of the inner bottom is 2 m and the maximum penetration does not sufficiently exceed to rupture the inner plate. It is believed that some discrepancies may also arise because the non-linear Coriolis and centripetal forces due to added inertia as well as the history of motions and associated radiation damping forces are neglected in the simplified FSI model. Such effects may lead to a slight increase in the time-period of the heave and pitch motions (see Figs. 14 and 17(a)) as demonstrated by LS-DYNA simulations. In any case, the grounding dynamics of centerline and off-center grounding events compares well with the hydrodynamic damping, restoring forces, damage extents, deformation energy and contact forces (see . A shortcoming of the proposed method is the accurate estimation of the lateral forces. This is based on the relative velocity between the ship and the rock used in the simplified FSI model (see Eq. (13)). LS-DYNA simulations demonstrate  that lateral force discrepancy may emerge from the longitudinal member of the hull bottom. In an off-center grounding of the box-shape ship, the longitudinal stiffeners close to the base of the rock experience transverse deformation (see supplementary videos of box shape ship and Fig. 19a). However, the present model of Simonsen (1997a) does not consider the transverse deformation of the longitudinal members; as a result, the ship approaches the centerline girder due to increased rolling and yawing motion as the longitudinal stiffeners are unable to provide sufficient transverse resistance (see the deformation of longitudinal girder in a resistance force plot in Fig. 19(b)). As a result, the path of damage is more straight in LS-DYNA simulation in comparison to the simplified model (see Fig. 19).

Numerical validation of a passenger ship centerline grounding
The bottom structure details and hull profile of the FLOODSTAND SHIP B (Luhman et al., 2019) vary along the length of the ship. Structural particulars are illustrated in Fig. 13. The results of the centerline grounding simulation are shown in Figs. 20-22. The total length of the damage estimated by the rapid FSI model was 82.6 m. This is 20% lower than the one predicted by LS-DYNA (see Fig. 20) while differences are also evident in heave and pitch motions that associate with the forces leading to structural deformation. For example, the vertical force predicted by the simplified method has significantly higher values during the first 9 s of the simulation (see Fig. 21). The inconsistency in vertical force arises from the distinct shape of the bow and fore shoulder regions of the ship. The differences (see Fig. 22), may be attributed to the fact that the simplified FSI model accounts for a simplified tearing mode for the dynamics of the flat bottom plate (Simonsen, 1997a). In reality, the bow and forward end regions are not flat. This results in a larger gap opening of the bow region and an overestimation of plate resistance and girder forces (see Appendices B and C). Vertical forces influenced by complex geometries are not accounted for by the simplified FSI approach. The mean deviation of the horizontal force is less than 0.5%. The average difference of vertical forces between the simplified FSI model and LS-DYNA simulation is approximately 28.5%. The total sliding energy from the rapid FSI model is marginally higher around 5.2%. On the contrary, the total deformation energy is almost 12.5% lower as compared to LS-DYNA results (see Fig. 21). Table 5 and Fig. 22, display a comparison between the damage extents predicted by LS-DYNA simulations, the simplified FSI model and the Simonsen approach for the case of centerline grounding. The longitudinal location (position) of maximum damage width obtained using the simplified model and LS-DYNA simulation are at 4.5 m and at 103 m from the bulbous bow tip, respectively (see encircled region in Fig. 22). This deviation is due to the complex bow geometry of the ship where the simplified model overpredicts the deformation width. This is not the case for the box shape ship (see Sections 4.2 and 4.3) where maximum and minimum damage widths obtained from the simplified FSI model and LS-DYNA simulations almost coincide.
The double bottom configuration along the length of the passenger ship varies. From the tip of the bulbous bow to a distance of 43 m towards the aft of the ship, the height of the double bottom is greater than 2.5 m. After 43 m from the bow tip, the double bottom height reduces to 1.6 m (see Fig. 13). The variations of resultant deformation forces for plates, girders, stiffeners, and floors are demonstrated in Fig. 22 and the supplementary material submitted for the case of passenger ship grounding. The simplified FSI model illustrates good quality results as it reasonably captures the forces  of the structural components. Moreover, as the position of the ship center of gravity reaches closer to the rock the pitch motion reaches zero (see Fig. 20). Once again, Simonsen's approach overestimates the damage length. As the damage length reaches approximately 145 m, the rock interacts with the girders (see Fig. 13). This increases the girder force as seen in Fig. 22c.

Numerical validation of a passenger ship off-center (B/4) grounding
Results for the case corresponding to off-center (B/4) grounding of a passenger vessel are shown in Figs. 23-26. The results of the surge, heave, and pitch motions match well. The total damage length of the ship when using the simplified FSI model is almost equal to LS-DYNA simulation. Comparison of the surge velocities from LS-DYNA simulation and present model shows that the grounding event terminates after 30.6 s and 33.2 s respectively. A constant surge velocity is attained after 30 s (see Fig. 23). The difference in the surge velocity may be attributed to the lower estimation of the horizontal force by the simplified FSI model. The simplified FSI model shows on average 10% lower values of longitudinal resistance in comparison to LS-DYNA (see Fig. 24). On the other hand, the trends of vertical force estimates compare well. The variation of vertical force is also reflected in heave motions where the simplified FSI model has slightly less heave amplitude.
The differences in the roll response between LS-DYNA and the simplified model are significant. LS-DYNA simulations show that when the rock is in between two longitudinal girders, roll motion is steadily increasing but is delayed in comparison to the values predicted by the simplified FSI model (see supplementary videos of passenger ship off-center grounding). During the first 6.7 s of the LS-DYNA simulation, the centerline longitudinal hull girder (Girder-1) is in contact with the rock (see Fig. 25). Transverse forces from Girder-1 subsequently incline the ship towards Girder-2, which is located nearly at the extreme width of the ship (see supplementary videos of passenger ship off-center grounding, Figs. 24 and 25). However, after 6.7 s the rock conforms to Girder-2 which bends. This ultimately affects ship sway and yaw motions (see Fig. 23). The trend of the sway directional lateral force predicted by the simplified FSI model is different from that of LS-DYNA. The differences in average forces predicted by both models are within the acceptable numerical error (just 2%). However, the differences in total deformation and sliding energy predicted by the simplified FSI and LS-DYNA models 10% and 4% respectively (see Fig. 24). The details of damage extents for the case of the passenger ship off-center grounding are shown in Table 6. The maximum deformation width and hull penetration estimated by the simplified FSI model are lower in comparison to those from LS-DYNA. Notwithstanding this, the differences in the maximum extent of damage is within 10% of LS-DYNA simulations. In both cases, the inner bottom ruptures (see Fig. 26). For off-center grounding when the ship rolls the location of the rock gets closer to the centerline and the longitudinal girders in between the rock tip and centerline girder may deform (see Fig. 26 after 94 m of the ship's length). Overall, the Simonsen technique underestimates total deformation energy and forces, resulting in an increase in damage length and the inability to approximate the real scenario. Table 7, summarizes the comparison of the computation time of the present simplified FSI model against LS-DYNA simulations. LS-DYNA simulations were performed using parallel execution with 120 parallel processing (MPP) processors. Table 7, shows that the simplified FSI model drastically reduces the computational effort in comparison to the fully coupled FSI model of LS-DYNA and can be used for stochastic assessment of ship damage extents.

Conclusions
This paper proposed a rapid FSI method for the realistic assessment of ship grounding dynamics. The model is validated by direct comparisons against LS-DYNA-MCOL and the simplified crashworthiness model of Simonsen (1997a) for the case of a barge and a modern cruise vessel. The new rapid assessment model provides more realistic results in comparison to the method of Simonsen (1997a) that assumes constant plate split angle and therefore overestimates ship damage extents. Comparisons against LS-DYNA-MCOL simulations show good agreement in terms of longitudinal and vertical     hull (e.g., hull curvatures towards ship bow). To resolve this discrepancy future research could focus on the development of an analytical formulation that accounts for transverse deformation of longitudinal members and tearing of the bow region of the ship. Notwithstanding this, the method may be considered significantly advanced in comparison to simplistic models and empirical data sets used by existing damage stability regulations. It also accounts for the influence of evasive maneuvers on the response (Conti et al., 2021). In this sense it could also serve as useful update in algorithms used by modern navigation decision support systems (Gil et al., 2020;Zhang et al., 2021Zhang et al., , 2022.

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.

Rigid body ship kinetics
According to Fossen (2011) and Lewandowski (2004) for the case of a symmetric displacement monohull the inertia and added inertia matrices are expressed as: where, m is the mass of the ship; I X , I Y and I Z are the moments of inertia of the ship around the x, y, and z-axes. In this paper the added mass coefficients in M A were obtained from linear strip theory. In the method proposed in this paper, the influence of time variation of added mass in way of contact was ignored and therefore the encounter frequency was assumed to be infinite during grounding contact. The mass multiplier terms in a rigid body inertia matrix (M ) arise because the origin is located at a distance ⃗ r G away from CoG. For the same reason, Coriolis and centripetal forces also include the effect of origin shift so that: Only small out-of-plane motions were considered. This means any changes in the waterplane area were independent of heave, roll, and pitch motions. The restoring force matrix given in Eq. (19) shows the combination of the effect of ship weight and buoyancy.
The damping ratios ζ w , ζ φ and ζ θ were evaluated numerically from the nonlinear time-domain 6-DoF potential flow solver of Matusiak (2021) following a numerical decay test. The validation of external ship dynamics is presented in Taimuri et al. (2020).
As part of the numerical decay test, the ship was assumed to progress in calm seas along a straight line with different initial velocities. She was released from an initial heave displacement, and initial roll and pitch inclination. The resulting heave, roll and pitch motions were measured separately to compute the damping ratio ζ , from the successive maxima (R n and R n+1 ) of the decaying motion given in Eq. (20).
Damping forces were then evaluated from the natural (heaving, rolling and pitch) frequency of the ship as presented in Eqs. (21) and (22). The curve fits were performed on the damping ratios (ζ w, ζ φ , and ζ θ ; Fig. 27) as a function of the Froude number. These values of damping ratios are used in a time-domain maneuvering simulation to evaluate radiation damping of the ship Eq. (21). The plots of decaying inclination (roll and pitch) and displacement (heave), and curve fitted damping ratios are illustrated in Fig. 27. The hydrodynamic damping coefficients B H were assumed coupled only for in-plane (surge, sway, and yaw) motions (Taimuri et al., 2019). Heave, pitch, and roll damping effects expressed by the term B p accounted for wave radiation effects.

Appendix B
Membrane straining of plate (Simonsen, 1997a) The membrane strain of the plate is illustrated in Fig. 4, which is defined as: where v 0 is half the length of the maximum transverse opening; u 0 is the longitudinal gap opening of the plate.; θ is the plate split angle; B de represents the deformation width and α is the flap angle representing the formation of the hinge line.
Description of the mechanics of transverse damage extents (Simonsen, 1997a) During ship grounding the stiffeners along the longitudinal direction of the bottom plate structure will deform. In addition, because of the influence of membrane stresses and transverse floor deflections deformation effects will extend to the nearest longitudinal. The critical distance of this longitudinal from the rock is defined as: where R R is the rounded tip radius of the cone representing the rock tip (Fig. 3); ϵ cr is the critical rupture strain of the material. Consequently, the penetration due to fracture is defined as: where, δ fra is the critical penetration causing the plate to rupture; s lf is the lateral extent of deformation from rock tip to the first longitudinal stiffener outside s l,c and ψ C is the plate wrapping angle on the rock.
The extent of the transverse deformation before fracture (intact condition) and the width of the damage in case of fracture were defined as: where, δ out is the penetration of the rock inside the hull. It is noted that the penetration into the hull δ out originates from the contact detection approach (see Section 3.4). (Sun et al., 2017) Deformation initiates once the base of the rock is in contact with the floor bottom plate (see Fig. 28). The vertical height of the deformed floor H def is defined as:

Description of the mechanics of longitudinal damage extents
where ∆ t is the deflection of the transverse web floor; φ is the cone semi apex angle; χ t is transverse floor spacing During deformation the floor is longitudinally transfers from line AE to line ACF Fig. 28. Therefore the critical fracture strain of the deformed floor is defined as: where b l,Gir is half the spacing in between two longitudinal girders; θ t is half of the arc CD Fig. 28, which is obtained from the geometry of the floor deformation as: The algorithm developed is illustrated in Fig. 29. The mechanics of intact and fractured plate conditions are specified in Eq. (25). The overall plate topology (i.e., unstiffened plate, girders, floors, and stiffeners) was assumed to remain in contact with the rock without fracture when the rock is adequately blunt (i.e., the penetration is lower than δ fra or the cone apex angle φ > φ c ). The critical apex angle φ c was defined according to Simonsen (1997a) as:

Plastic resistance forces of structural members
According to the upper bound theorem of plasticity, the resistance forces were computed by assuming a mode that accounts for the effects in way of the crack tip, bending hinge line, and membrane deformation. Appendix C summarizes the plastic intact and fracture forces defined by Simonsen (1997a) and Sun et al. (2017).

Appendix G
A comparison of the simplified contact model against LS-DYNA simulations is illustrated in Fig. 31. LS-DYNA considers the complete contact in between the surface of the hull and rock. This is illustrated in Fig. 31(a), where the hull bottom touches the generator of the cone at some distance from the tip of the rock. In the simplified FSI model the contact is detected once the tip is inside the hull (i.e., penetration into the hull becomes δ out ). After the first contact, the deformation mode follows the pattern described in Section 3.3 and Appendices B and C.

Appendix H. Supplementary data
Supplementary material related to this article can be found online at https://doi.org/10.1016/j.jfluidstructs.2022. 103589.