A Fully Coupled CFD-DMB Approach on the Ship Hydroelasticity of a Containership in Extreme Wave Conditions

: In this paper, we present a fully coupled computational ﬂuid dynamic (CFD) and discrete module beam (DMB) method for the numerical prediction of nonlinear hydroelastic responses of a ship advancing in regular and focused wave conditions. A two-way data communication scheme is applied between two solvers, whereby the external ﬂuid pressure exported from the CFD simulation is used to derive the structural responses in the DMB solver, and the structural deformations are fed back into the CFD solver to deform the mesh. We ﬁrst conduct a series of veriﬁcation and validation studies by using the present CFD–DMB method to investigate the global ship motion, vertical bending moments (VBMs), and green water phenomenon of the ship in different regular wave conditions. The numerical results agreed favourably with the CFD–FEA model and experimental measurements. Then, the extreme ship motions are studied in focused wave conditions to represent extreme sea conditions that a ship may experience in a real sea state. According to the conclusion drawn from the numerical simulations, it is founded that the focused wave case will lead to the increase of the longitudinal responses of the hull compared to regular wave condition, i.e., the heave, pitch, and total VBMs rise about 25%, 20% and 9%, respectively. In focused wave conditions, intensive ship responses and severe waves cause stronger slamming phenomena. It is found that the instantaneous impact pressure from the focused wave is higher and sharper compared to the regular waves and comes along with the obvious green-water-on-deck phenomena.


Introduction
Ships can experience slamming impacts while operating at high speeds, in severe sea states, or both. Severe slamming events can result in not only enormous local impact pressure but also hull girder global whipping responses, which is critical for the ultimate strength evaluation of a ship structure. El Moctar [1] presented a statistical report from the International Union of Marine Insurance [2] stating that approximately 36% of ship losses of all total losses, between 2001 to 2015, were associated with ships encountering harsh weather.
An increase in ship size tends to make the hull more "flexible" and leads to the wave encounter frequency being closer to the ship's wetted natural frequency. In such resonance frequency, comparatively small vibration-induced loads are superposed to the wave-induced load that may significantly violate the ship motions and deform the ship structure. This phenomenon is called the ship's hydroelastic effect, which is particularly important in determining design wave loads and structure strength evaluation of large ships [3]. The prediction of slamming loads on a deformable ship is a very challenging subject due to the strongly nonlinear free surface and the structure elastic behaviour. In such cases, the deformed ship structure interacts with the surrounding flow fields, which form a fully coupled fluid-structure interaction (FSI) system against the traditional rigid ship assumption. Hirdaris [4], Jiao [5], and Bakti [6] also supported this point of view, showing that the rigid body model may lead to inaccurate predictions of hydrodynamic loadings as well as the global ship motions. Therefore, as suggested by Hirdaris [4], a fully coupled FSI approach was used in this study to predict the flexible ship motions and hydrodynamic loads of an S175 containership in extreme wave conditions.
Over the past decade, the theoretical and numerical methods for predicting ship seakeeping performance based on hydroelasticity theory have gained momentum. To account for both fluid and structural features, the partitioned approach was adopted most commonly, in which the wave-structure system was divided into the fluid and structure parts and solved iteratively. The pioneering work of Bishop [7] developed a linear FSI model on the basis of 2D potential flow theory and a linear beam model on a monohull. Within this framework, the flexible structural characteristics of the ship was idealised as an elastic beam and interacted with the fluid forces, which were calculated from the strip theory. The theory was subsequently extended by Bishop [8] to a generalised beam model for floating vehicles or more complicated shapes. The three-dimensional hydroelasticity theory was established in the early 1980s by Bishop [8]; since then, a great deal of progress has been made in the development and applications of the hydroelasticity analysis on ship structures based on the potential theory [6,[9][10][11][12][13].
The potential flow theory provides an efficient insight into the initial design stage of shipbuilding; however, it cannot reproduce some important physical phenomena, such as wave breaking and viscous effects. As El Moctar [1] stated, the use of potential flow theory was questionable when applied to extreme wave conditions because the viscosity effects become critical for modelling large free-surface elevation waves. Therefore, the fully nonlinear computational fluid dynamics (CFD) method, which makes use of the Reynolds averaged Navier-Stokes (RANS) equations, is commonly used to investigate ship seakeeping in harsh seaways as an alternative. Extensive research has been proposed and validated the performance of the CFD method on the studies of ship seakeeping performance under different conditions-for example, coupled ship motions in deep water [14,15], ship resistance in restricted shallow water [16], and ship manoeuvring in irregular waves [17], to name a few.
The majority of current CFD studies on ship seakeeping were proposed by modelling fluid flow around a rigid ship. Although global ship motions and external loadings can be routinely obtained by CFD simulations, the hull section loads, e.g., vertical bending moment and shearing force, used for wave loads analysis cannot be directly obtained. Therefore, the CFD method can be further coupled with a structure solver to study the dynamic responses of a flexible ship on a free surface. In some cases, CFD and the finite element analysis (FEA) are coupled for the hydroelastic simulations based on a one-way coupling approach [18] or a two-way coupling approach [19]. The differences between one-way and two-way coupling methods on ship wave loads and hydroelastic responses were investigated by Lakshmynarayanana [18] and Takami [20] by using a commercial co-simulation interface (Star-CCM+ & ABAQUS). According to their results, the oneway coupling approach should not be used for the studies of the ringing and whipping phenomenon because the added mass effects in the elastic deformation are not accounted for, and it may result in an underestimate of the high vibratory components as well as the wave loads. Generally, there are three types of structure modelling of ship in the FEA: the pure-beam model, the beam-shell model, and the full ship model. First, the pure-beam model was built based on beam theory, i.e., Euler-Bernoulli beam or a Timoshenko beam, in which the ship's elastic deformations are computed by solving the beam equation of motions and external forces. El Moctar [1] proposed a study based on the Timoshenko beam model and coupled with a CFD model by a two-way coupling method to investigate the wave-induced structural loads of three containerships in regular and deterministic wave sequences. Their numerical results including the dynamics loads and structural vibrations agreed well with the experiments, which assessed the feasibility of using the beam model and transient RANS solver to obtain the short-term statistical measures of nonlinear ship responses. Secondly, the beam-shell model requires a more complex design in FEA, in which the ship geometrical models in the CFD and FEA solvers should be topologically equal. The wet hull surface is generated by using the shell elements and is kinematically constrained with the beam nodes in such a model. Prior to the FSI simulations, the structural modal analysis has to be applied to calibrate the beam profiles [5,19]. Such an FEA model was coupled with a CFD solver based on a two-way coupling approach through co-simulation interface (Star-CCM+ & ABAQUS) by Lakshmynarayanana [21], Lakshmynarayanana [22], and Jiao [5] to predict the motions, wave loads, and hydroelastic vibrations of an S175 containership in waves with forward speeds. Their results, including global ship motions and vertical bending moments in different regular wave conditions, were comprehensively analysed and validated with the existing experimental results [23]. The third method is the full ship model, which discretises the whole ship and substructure by using the finite shell elements in the FEA software. This method is capable of accounting for the complex deformation behaviours of the ship and predicting local stress distributions; however, it is computationally very demanding. Ma [24] modelled the full ship with sandwich composite material by using ANSYS and coupled it with CFX through the ANSYS workbench. Their FSI results were further studied with stress analysis, which revealed that the most vulnerable region in the sandwich plate structure was the core-skin interface near the girder.
However, the abovementioned CFD-FEA coupling method requires significant computational efforts [6,21], making it difficult to extend applications to more complex wave conditions. Recently, a CFD-DMB method was first proposed by Lu [25] on the investigation of the hydroelastic behaviour of a large VLF. The DMB model has a similar setup with the traditional hydroelastic experiments. In this approach, the ship hull is divided into multiple floating rigid sections, whereas a stiffness matrix based on an Euler-Bernoulli (EB) beam theory, equivalently representing the structure's stiffness, is employed to connect ship neighbouring sections. In this method, only a one-dimensional beam is modelled in the structural solver, and it can significantly save computational time and effort compared with the CFD-FEA method. Based on the above design, the dynamic motion of the ship is both affected by the hydrodynamic forces from fluid solver and restricted by the deformation conditions of the equivalent beam properties from the structural solver. Such a CFD-DMB coupling approach has been applied successfully in the literature-for example, flexible wind turbine blades [26] and aircraft wings [27]. However, in these applications, a single-phase solver is applied to calculate the air forces at the geometry surface with a low mass ratio. For multiphase cases, the DMB method has been successfully employed by Li [28] on the dynamic motion of a close-loop WEC array in waves. Lu [25] investigated the hydroelasticity of very large floating structures (VLFs) in waves, which was solved in both frequency [29] and time domains [30]. Recently, Bakti [6] studied the forward speed effects on the elastic responses of an analytical Wigley hull by using a BEM-DMB coupling method in both regular and irregular waves. According to their results, a noticeable increase of ship vertical bending moments occurs when the ship speed increases due to the whipping phenomenon.
Most ship hydroelastic research in the literature was presented in regular wave conditions. However, regular waves are unrepresentative of actual extreme events, as stated by Ning [31]. On the other hand, random wave simulations require very long runs in order to capture the near-extreme events because they are so rare in any random wave series. In practice, the focused wave theory is used as an alternative to a long-term irregular wave. The concept of the focused wave was first proposed by Davis [32] by modulation of a series of conventional sinusoidal wave trains generated from a prescribed wave spectrum and superimposing the crests. Since then, experiments and numerical investigations have been carried out by using the focused wave groups to replace the irregular waves in the nonlinear wave-wave interaction studies. For example, T. E. Baldock [33] created wave focus events through the superposition of regular wave trains based on the linear wave theory and investigated the effects of nonlinearity of wave-wave interactions. Ning [31] studied the propagation of the focused wave groups with different incident wave parameters and compared the solutions between the first-and second-order wave theory. Recently, the focused wave theory has been adopted to investigate wave-structure problems and predict extreme loads on offshore structures. These studies were initially focused on a simple floater [34], and then extended to complex floating structures, such as wave energy converters (WECs) [35] or semi-submersible wind turbine foundations [36,37]. Zhou [38] and Zhou [39] further extended their applications of studying the dynamic motions of a fully coupled aero-hydro-moored version of a floating offshore wind turbine (FOWT) subject to focused wave conditions.
To the best of the authors' knowledge, the research using transient FSI simulations to investigate ship hydroelasticity in extreme wave conditions is very limited in number and scope. The reasons for this mainly stem from the computational burdens of determining maximum values of ship responses in irregular waves through direct simulation by using Navier-Stokes equations. In this study, we overcome such issues by studying ship hydroelastic behaviour and slamming loads in extreme waves based on the CFD-DMB method and focused wave theory. It is expected that the results obtained from this improved and validated numerical tools can provide a more precise and more detailed insight into the physical phenomena of the ship dynamic motions and its hydrodynamic loads in real sea states. The results proposed in this study could also help to assess the structural integrity of the ship longitudinal strength, which serves as an improved technique by which to evaluate unconventional ship designs.
The rest of this paper is organised as follows. In Section 2, the numerical setup and the methodology used in the present study are discussed. In Section 3, the verification and sensitivity studies are presented for focused waves and the CFD-DMB model. In Section 4, the numerical results for the study on the dynamic motions of a flexible S175 ship in regular wave conditions are presented. In Section 5, a comprehensive comparison of the ship global motion and hydroelastic behaviours between focused and regular waves are presented. The discussion and conclusions are drawn in the final section.

Problem Statements
In this study, the numerical simulations were performed for an S175 type containership with a scale ratio of 1:40, as displayed in Figure 1. This type of ship geometry has been commonly used in the literature for ship hydroelasticity research. These studies include the BEM-FEM coupling research by Kim [12], the CFD-FEA coupling research from [5,[18][19][20]40], as well as some experimental results as reported by Chen [23]. Only the bare hull was modelled for the seakeeping investigations in this study. The main dimensions and the characteristics of the full and model scales are listed in Table 1.  In this research, an array of six pressure gauges (P1, P2, P3, P4, P7, P9) were placed on the bow flare and bow bottom areas for bow wave pressure measurement to investigate the impact wave loads on the ship when sailing in harsh weather conditions, as shown in Figure 2a. Moreover, an array of three pressure gauges (P5, P6, P8) were arranged at the side of the hull to measure the slamming pressure at the starboard as shown in Figure 2b.

Numerical Modelling
In the present FSI study, a two-way CFD-DMB coupled method was implemented by coupling an open-source CFD toolbox OpenFOAM-v2012 [41] and a structural solver MBDyn-v1.73 [42] with an in-house data communication tool. The main features of the fluid, solid solver, and coupling method are illustrated in this section.

Fluid Dynamics
The mass continuity and momentum equations for a transient, incompressible, and viscous fluid can be written as ∇ · U = 0 (1) where U refers to the velocity of flow field, ρ is the mixed density of water and air, g is the gravity acceleration, P d refers to the dynamic pressure, τ is the dynamic viscosity, and F is the surface tension. The VOF technique [43] is adopted to simulate the free surface in the numerical domain by solving the following transport equation for the scalar quantity, a, which represents the volume fraction of fluid for each cell. In multiphase flows, the volume fraction is assigned a value of 0 when the cells are filled with air, and when it is 1 the cell is filled with water. The air-water interface is capable of capturing at the free surface cells with a volume fraction between 0 and 1: To reduce the numerical smearing and sustain a sharp interface between water and air, OpenFoam in particular adds an artificial compression term U r (1 − a)a (the third term in Equation (3)) [44], where U r is the artificial compressive velocity which only functions near the free surface due to the inclusion of (1 − a)a.
The multiphase solver interFoam used in this study is based on a finite volume method that solves the unsteady RANS equation in an iterative manner. The density ρ and the kinematic viscosity v of the fluid in simulation are reported in Table 2. In the numerical model, the convection terms were discretised with a Gaussian linear upwind scheme, whereas a second-order Gaussian linear corrected scheme was adopted for the diffusive terms. The turbulence model selected for the simulation is the shear stress transport (SST) model [45], which blends the best features of the near-wall accuracy of the k − ω model and the free-stream accuracy of the k − ε model. The temporal discretisation was performed with a first-order backward Euler scheme. PIMPLE (a combination of pressure implicit with splitting of operator (PISO) and the semi-implicit method for pressure-linked equations (SIMPLE)) was utilised to solve the pressure velocity coupling equations.

. Computational Domain and Boundary Conditions
The finite volume mesh was generated by using the OpenFOAM default mesh generation tool "SnappyHexMesh" based on cell-splitting and mesh-fitting techniques [41]. A uniform background mesh was initially generated and used to project and snap cells onto the geometry, and then the mesh refinements can flexibly be specified on edges, surfaces, and volumes to obtain optimum geometry feature resolutions. The numerical domain used in this study simulates ship motions in deep-water conditions, which extends in the three dimensions, i.e., −1.5 L < x < 2 L, −0.6 L < y < 0.6 L and −1.5 L < z < 0.5 L, where L refers to the ship length between perpendiculars. The grid density at the free surface is progressively refined until it fulfilled the guidelines from [46], in which a minimum of 100 cells per wavelength and 12 cells per wave height were used on the free surface modelling in this study, as shown in Figure 3a. To ensure that the high Reynolds number flow features are approximately captured, the grid density at the area around the ship hull is further refined several times with boundary layers, primarily maintaining the adjacent wall layer-thickness coordinate y+ close to 30. It is worth noting that the proper wall functions were implemented for the hull patch surfaces to model the approximate wall behaviour when the adjacent layer thickness stayed in a log-law region. The numerical domain with a wave height of 0.12 m and encounter wave frequency of 5.581 Hz is chosen as a representative case, which is displayed in Figure 3 with a general view of the computational mesh and mesh-refinement zones. The boundary conditions in the present CFD study are shown in Figure 4 as follows.
At the left boundary inlet, the velocity inlet was prescribed as the incident wave and current, whereas the pressure was set as a zero gradient. At the right boundary outlet, the current velocity outlet was applied to preserve the conservative of flux inside the computational domain. The boundary condition of the domain top part was set as the atmosphere. The lateral sides were set as symmetry planes to avoid wave reflection at the boundaries. The bottom boundary was set as a symmetry plane for deep-water modelling. The moving wall boundary condition with zero pressure gradient was defined on the surface of the ship hull.

Wave Generation and Absorption
An open-source toolbox "waves2Foam" was used in this study to generate and absorb free surface waves in the computational domain [47]. A relaxation zone technique was adopted to provide better wave quality and to remove spurious reflection from numerical simulations. Two relaxation zones were placed at the inlet and outlet boundaries of the numerical wave tank. Equations (4) and (5) specify the main functions of the relaxation zone technique. We have where φ R refers to either the velocity or volume fraction of water a, the superscript c and t defined the value of computed and target, respectively. The weighting function a R is always equals to 1 at the interface between the non-relaxed computational domain and the relaxation zones, χ R is a value between 0 and 1. The relations between χ R and a R are shown in Figure 5. Figure 5 shows how a series of surface elevation gauges (P1, P2 and P3) were placed inside the numerical domain to monitor the quality of wave generation.

Focused Wave Modelling
The focused wave group implemented in this study is based on the NewWave theory [31], which generates the extreme wave event from a specified sea spectrum by superimposing several relatively small amplitude waves at a chosen point and time. The spectral shapes of irregular waves were implemented by the JONSWAP spectrum [48]. The significant wave heights H s , peak angular frequency f p , and peak lifting factor γ are the main parameters of the JONSWAP spectrum. We have where, β is the coefficient related to the peak lifting factor γ, σ is the shape factor with a value of 0.09. For the linear NewWave theory, the amplitude of each wave component a i of frequency f i is defined as where S( f i ) is the surface spectral density, ∆ f is the frequency step (which depends on the number of wave components N and bandwidth), and A 0 is the target theoretical linear wave amplitude of the focused wave [31]. The extreme wave represented by linear NewWave theory is simply the scaled auto-correlation function corresponding to a specified spectrum. The linear free surface elevation η (1) and horizontal and vertical velocities u (1) and ω (1) are given by where, z is the free surface amplitude measured from the mean water level (MWL), x 0 , t 0 are the predefined focal location and focal time, respectively, a i is the amplitude of wave components, g is the gravity acceleration, h is the water depth, ϕ is the phase angle, is the wave number, and ω i = 2π f i is the frequency. Modulation of phase angle among individual wave components (Equations (9)-(11)) can superpose the wave peak at a fixed time and position with the mathematical representation as In CFD modelling, a uniform current is often applied at the boundary inlet to simulate the ship forward speed. However, the current effects influence the focused wave modelling, including wave steepness, focusing time and position, owing to the nonlinear interaction between the wave and current [49]. In the present study, the nonlinear interaction between the wave and current was not considered; therefore, a linear focus wave with the current model proposed by [50] was implemented as given in Equation (13). We have where V is the uniform current speed. Based on this wave-current dispersion relation, the circular wave frequency ω i and wave number k i can be calculated by using a mathematical iterative method, the positive root of Equation (13) will find the value of k i . The rest of the wave parameters, i.e., wavelength and frequency, are updated accordingly. The current at the boundary inlet provides a horizontal flow velocity and modifies the focused wave velocity from Equations (10)- (14). To achieve the peak of wave trains energy at a target time and position, the focal location and time are required to be modulated accordingly.
We have

Structural Dynamics
To model the nonlinear hydroelastic behaviour of the containership in question, the discrete module beam (DMB) method was applied by adopting a multi-rigid body approach in the structure solver MBDyn [42]. The concept of the DMB approach is similar to the standard hydroelasticity experiment's setup of [6].
The ship hull was divided into m sections, connected by m − 1 beam elements, consisting of a multi-body structure system. For each body of the system, the Newton-Euler equations of motion were established in the differential algebraic form as a set of firstorder equations together with the constraint equation, resulting in a system of differential algebraic equations (DAE), as follows. We have where M denotes the inertia matrix of the rigid body, x denotes the translational and rotational parameters in the global reference frame. P refers to the momentum of the body, λ is the vector of the Lagrange multipliers for the constraints, f is the external force and moment vector exerted upon the body which might be related to its displacement and velocity as well as time, φ is a set of kinematic constraints applied on the body, and φ T x is the Jacobian of ϕ with respect to the generalized coordinate.
The distance between the centre of gravity of two adjacent sections is taken as a spatial beam for considering the effects of structural deformation, as shown in Figure 6. The beam element used in the present study is a three-node beam element implemented in the MBDyn software by a finite volume approach for the multibody formulation of three-node beam elements based on the geometrically exact beam theory (GEBT) [51]. The internal forces and moments were evaluated at the cross-section at the evaluation points (shown in Figure 6) and related to the geometrical strains and curvatures via the constitutive law with the following equations. We have where F xx is the axial force component, F yy and F zz are the shear force component, M xx is the torsional moment component, M yy and M zz denote bending moment components, ε xx , γ yy and γ zz are axial strain and shear strain coefficients, κ xx , κ yy and κ zz are the bending curvature parameters, and f denotes as an arbitrary function of beam material constitutive law. The performance of the above beam model has been validated by Bauchau [52], Liu [53] and Liu [26]. According to Bauchau [52], this beam model is valid when the structure's plane section remains plane (no shear deformation) and deformed beam angles are small. The effects of cross-section warping are assumed to be small and neglected.

Ship Modelling
The deformable ship hull was divided into 20 sections whereas a stiffness matrix based on the Euler-Bernoulli beam theory, equivalently representing the structure's stiffness, was employed to connect the neighbouring sections. In its numerical representation in MBDyn, the distributed mass of each ship section was modelled as a lumped mass point located at its centre of gravity according to the mass distribution data from [19], as shown in Figures 7  and 8. A uniform beam stiffness was applied with vertical bending stiffness, which is shown in Figure 8. The mass moment of inertia, I * yy , I * zz , and I * pp of each ship section was calculated based on the simplified approach from Bakti [6]. In their approach, the moment of inertia of each section (I y ,I z ,I p ) was calculated by the product of the first moment of area and the centroidal distance of the area from the given axis. The parallel axis theorem was applied when the centre of the ship cross-section was not coincident with the beam node. Afterward, the moment of inertia can be calculated by using geometric properties of the cross-section by multiplying a factor of a, such that I * yy = a * I y . The factor a = 8.0 was used in this study, which was determined by matching the first two natural frequencies of the vertical bending modes of the beam with the values calculated by FEM. The material property of the beam was defined as steel, i.e., elastic stiffness E = 210 Gpa and Poisson ratio v = 0.3. In this study, the structural damping was estimated approximately as 1% of the critical damping, as recommended by [8].
The DMB model in the present study was constrained for the y-axis (no sway) translation and rotation about the x-axis and z-axis (no roll and yaw) by imposing total joint elements among beam nodes. Moreover, another set of total joint elements was employed to restrict the DMB model in the direction of the x-axis to prevent the longitudinal drift caused by the wave and current. Based on the above design, the elastic deformation of the ship hull is affected by the equivalent beam properties and restricted by the constraint equations.

Modal Analysis
Before the hydroelastic computation, the modal analysis of the DMB beam model was conducted by using arpack solver in MBDyn software to provide information about beam vibration behaviour, such as the natural frequencies and modal shapes in the dry condition. The beam was assumed in a vacuum condition, and the elastic behaviours were considered by the multibody dynamic beam theory with equivalent beam stiffness. The beam section profiles (including section dimensions and thickness) were calibrated by matching the first three modes' natural frequencies with the natural frequencies of the real ship from experiment [5]. The obtained natural frequencies for the first three orders of the vertical bending mode of the model in the dry condition are summarised in Table 3, and the corresponding mode shapes are shown in Figure 9. After the calibration, the closed beam profile with a closed cross-section of 0.08 m × 0.05 m × 0.005 m was used throughout all later cases in this study.

CFD-DMB Coupling Approach
To approach a fully coupled CFD-DMB simulation, there are two key components to be considered: (1) coupling methods for data communications and data mapping, (2) dynamic mesh handling to update structural motions in CFD mesh. In this section, the present CFD-DMB coupling framework will be discussed based on the above two aspects.

Coupling Algorithm
The two-way FSI algorithm was applied for data communications between the fluid and structural solver to satisfy the geometrical compatibility and the equilibrium conditions on the interface. Generally, a two-way algorithm can be further divided into strong and weak coupling methods. The two-way strong coupling scheme, indicated by Benra [54], has the second order's convergence rate, which is appropriate for the FSI problems with a strong dependence between fluid and structure. Such a two-way strong coupling method has been applied for the flexible ship seakeeping problem by using the commercial co-simulation interface (Star-CCM+ & ABAQUS) or open-source coupling library PreCICE [55].
As shown in Figure 10a, the two-way coupling method has a complex workflow, including PIMPLE and outer iterations in each time step. The convergence of the PIMPLE iteration has to be fulfilled first, and then the data communication between the fluid and structure solver will be executed until the outer iterations converge. The strong coupling method still poses difficulties in the present stage due to the large demand for computational effort. Therefore, this study used the weak coupling method as an alternative because it has a simplified workflow, as shown in Figure 10b. As can be seen from the figure that there are no outer iterations required inside the time step; therefore, the communications between the fluid and structure solver can be performed efficiently, which results in a significant reduction in the overall computational cost. A detailed workflow of the present CFD-DMB approach is shown in Figure 11. At the start of the time step, the PIMPLE iteration is executed first to solve the RANs equations of fluid fields. Until the PIMPLE iteration reaches convergence, the calculated fluid forces and moments at patch surfaces are transposed to the corresponding beam nodes in the MBDyn solver. Then, MBDyn solves the differential equations (Equations (15)- (17)) to compute the structural responses and the structural deformations are fed back into the CFD solver to deform the mesh.

Dynamic Mesh Motion
To update the fluid mesh motion according to the solution from the structure solver, the dynamic mesh-handling technique of mesh motion was applied in this study based on the finite volume-based Laplacian mesh motion for the displacement field. This type of dynamic mesh operation solely involves the displacements of the mesh points without altering the topological information of the mesh [56]. The distortion of internal parts of the solution domain relies on the solution of a diffusion (Laplace) transport equation for the displacement point fields using the equation shown below, where γ is the displacement diffusion coefficient, and d is the point displacement field. The rate of diffusion of the displacement from the boundaries to the internal mesh regions is defined with the spatially varying diffusion coefficient y given as a function of the distance between the point and the mesh boundary γ = γ(r). In order to calculate the displacements at the mesh points (cell corners), interpolations from the cell-centred values to mesh points are required. The mesh motion boundary conditions are applied at the moving bodies, using data computed based upon the calculations from the structure solver MBDyn.

Verifications
The following section will present two verification studies. In Section 3.1, a sensitivity study is presented on the focused wave generation based on the NewWave theory (Section 2.3), and the numerical results are compared with the analytical solutions. In Section 3.2, a verification study is carried out on the effects of mesh grids and time-step sizes on the present CFD-DMB coupling method at the wave-ship resonance condition (λ/L = 1.2, H = 0.12 m). The results, including dynamic ship motions (heave and pitch) and vertical structural loads (VBMs), are analysed in order to justify the correctness and accuracy of the present FSI model.

Sensitivity Study of Focused Wave Generation
The presence of uniform flow poses additional difficulties in modelling-focused wave groups in CFD due to the wave-current interaction. Similar numerical studies presented by Zhang [49], Markus [57] and Li [50] proved that currents significantly change the focused wave elevations, peak wave period, focal time, and position. However, such wave-current coupling effects are not one of the aims of this study; the focused wave group with current was generated based on the modified NewWave theory discussed in Section 2.3, and the focused wave height and peak wave period were then calibrated accordingly to meet the desired focused wave shape.
In this verification study, a 2D numerical wave domain extends into two dimensions of −9 m < x < 11 m, −6 m < z < 2 m. A focused wave group was generated with the significant wave height of H s = 0.12 m, and the wave peak period is chosen as T p = 1.78 s in the model-scale. The uniform current speed was applied at the boundary inlet with a value of 1.80 m/s, which implies the ship's forward speed of Fn = 0.275. The focal position was set as x = 0.0 m, and the focal time was set to be 15.0 s. The frequency band between 0.125 Hz to 2.0 Hz and 50 individual wave components were used to generate this focused wave group. The mesh information of the mesh convergence study was summarized in Table 4. To verify the quality of focused wave modelling, the discretisation error from mesh grids and time steps are studied and summarised in Tables 5 and 6. The numerical results of wave elevations are further compared with the analytical solution of a focused wave group without current as shown in Figure 12. The figure shows that the simulated focused wave preserves favourable wave shapes, and the peak wave elevation is focused at the focal position and time, as desired. With the consideration on the computational cost, it was concluded that the medium size of mesh and time step of 0.0015 s was the optimum choice to generate the waves in this study.

Sensitivity Study on Flexible Ship in Regular Waves
Prior to present analysis of the CFD-DMB model of ship responses in different wave conditions, it is necessary to conduct related convergence and uncertainty analyses of the FSI model to verify the coupling framework performance. As Huang [58] pointed out, the uncertainties in the modelling of fluid dynamics by a CFD solver are generally much greater than the uncertainties associated with the structural responses by the structural solver. Therefore, in this section, a verification study is conducted in the CFD solver by changing three different mesh densities and three time steps (details shown in Table 7).
The ship-wave resonant frequency condition (λ/L = 1.2, H = 0.12 m) was studied because large motions and loads tend to cause the largest numerical errors [15]. It is worth mentioning that the change of grid resolution was only applied on the free surface region, whereas the mesh discretization was not altered in the background mesh. The simulated time series of the heave and pitch motions and the VBMs at the amidship section (L/2 from FP) by the three different grid densities and time steps are comparatively shown in Figures 13-15. The time series of the flexible ship heave motions were monitored on the real-time displacement of the beam node on the ship at LCG. It is worth noting that the structure solver MBDyn does not have an internal rotational degrees of freedom [42]. The pitch angle is calculated by using the relative difference ratio between two nodes, at ship bow and LCG, based on a linear approximation equation, as given in Equation (20), where θ is the pitch angle in degree, subscript c points to the beam node at ship LCG, and b denotes the front node of beam. The raw pitch data is given in radians; therefore, a factor of 180/π is multiplied to convert it into degrees. Figures 13 and 14 show the time series of the heave and pitch signals from three different mesh densities and time steps, respectively. The results by different grid densities are very close, although the medium and coarse meshes slightly underestimate the heave peaks. Figure 15 shows the simulated time histories of the ship VBMs among different mesh resolutions and time steps. The data of wave-induced VBMs were filtered by the static still-water bending moment (SWBMs) from the total VBMs. The VBM signals in Figure 15 have shown similar shapes with the numerical results from [59], in which a strong asymmetry between sagging and hogging values exist due to the nonlinearity. To better understand the influence of waves on the high-frequency vibration characteristics, a fast Fourier transform algorithm (FFT) was applied to process the VBMs signals, and the results are shown in Figure 16. The main peak frequencies can be clearly observed from both figures, which implies the excitations from the encounter wave frequency, while the high frequency components can be captured up to the ninth order, which is denoted as the springing effects. It has to be noted that the high-frequency components are better recorded in fine mesh condition, but in general the extreme VBMs can be well estimated by using the medium mesh with less than 4.6% of the difference, as compared to fine mesh.
All simulations were performed by using the Archie-West UK (HPC) facility with 40 Intel Xeon 2.0 GHz cores. The average time for the regular wave case to obtain 20 stable periods consumed approximately 120 physical hours, and it will increase while running for focused wave cases.

Flexible Ship Responses in Regular Wave Conditions
In this section, the predicted global ship motions, i.e., heave, pitch and VBMs of the S175 ship model, are presented in different wave conditions based on the CFD-DMB method. The presentation of these results is fundamental for the subsequent focused wave analysis. There are five different incident wavelengths ranging from (λ/L = 0.9−1.6) with a constant wave height of H = 0.12 m, involved in this study. The detailed wave parameters are summarised in Table 8.

Monitored Incident Waves
The encounter waves for the ship sailing in sea state 6 with different wavelengths (i.e., λ/L = 0.9, 1.0, 1.2, 1.3, 1.6) measured at P2 are shown in Figure 17. As seen from the wave results, the time series of the monitored waves preserved shapes well with less than 6.1% of wave height dissipation while propagating, which confirms the capability of the developed CFD wave generation technique.

Ship Global Motions
The original time series of the flexible ship motions including heave and pitch among different wavelengths are shown in Figures 18 and 19. As can be seen in the figure, the vertical motion signals of the ship are all sinusoidal at any wave conditions that present the same sinusoidal characteristics of the induced waves. The largest heave and pitch motions are observed at the wave condition of λ/L = 1.4.  The heave and pitch motions are further calculated as response amplitude operators (RAOs) and plotted against the nondimensional parameter of wave/ship length ratio λ/L in Figures 20 and 21. To validate the numerical results, these RAOs are compared among different methods, including the numerical CFD-FEA models by Jiao [19] and Lakshmynarayanana [22], as well as experiments of Chen [23]. Figure 20 shows that the heave RAOs from the present CFD-DMB method generally show good agreement with the experiments by Chen [23] at short waves; however, the peak value at ship-wave matching resonance is captured at λ/L = 1.4 rather than λ/L = 1.2 as in the experiments. The correspondence between the present model and the CFD-FEA model is favourable in heave predictions except in the nondimensional frequency range of 1.4-1.6, where the magnitude of heave is overpredicted by about 25%. Figure 21 presents the pitch RAOs of the present method and the CFD-FEA model by Jiao [19], which produces small differences in the high-frequency region; however, there is an overestimation in long waves by an average of about 8.6%. Some differences in the pitch RAOs can be noted between the present model and the CFD-FEA model by Lakshmynarayanana [22] in long waves, especially with a peak difference at λ/L = 1.5 of approximately 15%. All numerical models overestimate the trim angle at long waves by approximately 25% compared with the experiments.

Wave-Induced Vertical Bending Moments
The longitudinal distributions of the total VBMs at each ship section estimated by the present CFD-DMB method are given in Figure 22 for different wavelengths. The VBM values are nondimensionalised by M/ρgL 2 Bξ, where ξ is the monitored wave amplitude. As can be seen from the results, the distributed VBMs show similar trends for different wavelength cases, and the peak crest and trough VBM values are found at the waveresonance condition (λ/L = 1.2), which implies most enormous bending moments occur at the wave-resonance condition. Both hogging and sagging VBMs show strong asymmetry behaviour along the ship longitudinal, the peak sagging VBM occurs at ship section 10, and the peak hogging VBM occurs close to ship section 12. The distributed VBMs at resonance condition (λ/L = 1.2) are selected as representative cases to further compare with the numerical CFD-FEA results from [5] as shown in Figure 23. The results revealed that the VBMs predicted from the present method show a similar trend to the results from the CFD-FEA method, especially the excellent agreement of the total sagging moments. However, the peak crest value of the hogging VBMs is detected at Section 10 in the CFD-DMB method, rather than at Section 11 in the CFD-FEA method. Figure 23. Comparison of the distributed VBMs at each ship section obtained by using the CFD-FEA method with those from [19]. Adapted with permission from [19], 2021, Jialong Jiao.
The VBM RAOS at the amidship is further compared among different numerical methods and experimental measurements. As shown in Figure 24, there is an underestimation of the results from the present CFD-FEA model at some wave conditions, i.e., λ/L = 1.0, 1.4 by a maximum of about 15%, compared to the results from Lakshmynarayanana [22]. These deficiencies may be caused by the implementation of a structural damping ratio, i.e., 0 in the CFD-FEA model and 0.01 in the CFD-DMB model. The damping effects on structural loads are studied by Lakshmynarayanana [22], which point out that the increase of structural damping from 0 to 0.01 may reduce the bending loads by about 25%. However, in the present study, applying structural damping to stabilize the structural solver and represent the material's real damping behaviour is necessary. Overall, the results give a fair depiction of the structural loads of the flexible ship to different wavelengths, enabling a further facile step for the focused wave study.

Green Water Visualisation
To better understand the green water phenomenon under different wave conditions, the physical views of the violated free-surface flow during slamming events at λ/L = 0.9 and λ/L = 1.2 are shown in Figures 25 and 26, respectively. Both figures include two phases of the ship motion state with four time steps within one wave cycle, from ¼ T to T. A severe slamming and green water on deck can be observed when the ship is sailing in λ/L = 1.2; however, such events are not noticeable at λ/L = 0.9 because the vertical motion of the ship is relatively low.

Comparison of Ship Motions and Impact Wave Loads in Focused and Regular Waves
In this section, the wave elevations, global ship motions and impact wave loads of the S175 model in regular waves and focused wave groups are presented in a comparative manner. In order to control the variables, the wave parameters, including significant wave height and trough-to-trough period of the focused and regular waves are set as identical. Moreover, to better understand the slamming impacts and green water on deck at extreme wave conditions, the flow field around the ship is visualised with the impact pressure predictions. Figure 27 compares the time-series free surface elevations of the focused wave group and regular waves (λ/L = 1.6) for the same significant wave height (H s = 0.12 m), troughto-trough period (T p = 1.4 s), and uniform current speed of U = 1.80 m/s. The focal location was predefined at the ship's LCG, and the focal time was set to be 15 s. These values were determined by considering the computational cost and arranging a sufficient time for wave propagation. The time records of 11.5 s to 19.5 s duration of the wave propagation were presented for comparison.

Wave Elevations Analysis
As the wave elevations are shown in Figure 27, the focused wave shows good agreement with the theoretical solution; however, a slight front shift of the focal time can be observed. The trough-to-trough period of the simulated focused wave group and regular wave is almost similar. As can be seen from the figure, the focused wave group tends to superimpose a higher wave elevation at the focal time, raising the wave elevation by about 25% from 0.061 m (regular wave) to 0.083 m. It is known that the increased wave elevation contains higher energy, which implies that forces due to focused waves are expected to be larger than that in the regular wave condition. The enlarged fluid forces further violate the ship response and corresponding elastic behaviour, which will be discussed in Section 5.2.

Flexible Ship Motions
The time histories of the heave motions induced by the focused wave are compared with that of the regular wave (λ/L = 1.6) as shown in Figure 28. As seen from the figure, the ship in the focused wave induced a significant motion response of the value of 0.095 m at t = 17.5 s, which increases the heave amplitude by about 25% compared with that in the regular wave. After the focal time, the heave motion decays progressively. The asymmetric behaviour can be noticed in heave signals, in which the trough of heave motion is generally much greater than the crest value, as the ship experienced obvious sinkage due to the forward speed and dynamic effects.  Figure 29 shows a comparison of the time series of the pitch motions between the focused and regular wave conditions. The pitch motion in the focused wave group shows a similar trend as the heave motion but with a phase shift, which reaches its peak after the focal time and decays progressively. As can be seen, the pitch amplitude increases by about 20% at the physical time of 18 s compared to that in the regular waves. The enlarged trim may result in severe whipping effects, which may endanger hull girder integrity. The time histories of the total VBMs at the ship amidship in the focused wave are compared to those in the regular waves (λ/L = 1.6), as shown in Figure 30. The sagging moment reaches its lowest value at the time instant of 16.5 s, with an increase of about 25% compared to regular waves. The total VBMs signal is further filtered by using a low/high pass filter to extract the wave-induced VBMs (green dot) and high-frequency VBMs (green dash dot). It can be seen from the figure that the high-frequency components were more pronounced between 16 s to 18 s, which is right after the focal time, accounting for about 20% of the total VBMs. The high-frequency components are mainly attributed to the nonlinear components in focused wave group and structures. To better assess the VBMs of each section along the ship length, the longitudinal distribution of VBMs of focused wave and regular wave are displayed in Figure 31. It is seen that the sagging VBMs is generally much larger than the hogging VBMs for both cases. The largest sagging VBMs appears at Section 10 in front of the largest hogging VBMs due to the contribution of nonlinear whipping loads components, which is evidenced by Jiao [19]. There is a noticeable increase by about 9.8% of sagging moments at Section 10 in the focused wave compared to that in the regular waves. The largest hogging moment is shifted from Section 12 (regular wave) to Section 11 (focused wave).

Slamming Loads and Green Water on Deck
In this section, the slamming loads and green water on deck phenomenon are analysed between regular and focused wave conditions. It should be noted that the time series of impact pressure data have a sampling frequency of about 667 Hz associated with the selected time step of 0.0015 s, which may not be sufficiently small to capture the instantaneous impact peak [5]; however, this sampling frequency is determined as a compromise between accuracy and the computational costs.
The time series of the impact pressure at front bow measurement points P1-P3 are compared between the regular and focused wave as plotted in Figure 32. The measurement points P1, P2, and P3, located at the centre line at the ship bow, capture the slamming impact loads and bow entry water pressure. The figure shows that the pressure signal at P1 generally has the sharpest peak due to its relatively small dead-rise angle compared with P2 and P3. However, the pressure signals contain high-frequency noises, which may be related to the ship ringing phenomenon. This part of the analysis will be carried out for future studies. The slamming pressure becomes higher and sharper in the focused wave at 15 s to 17 s compared to the regular waves, which indicates the ship experienced a severe slamming phenomenon.
The bow bottom gauges at P4, P7, and P9 are immersed in the water, which the hydrostatic pressure dominates. The pressure signals show general regular shapes with the occurrence frequency of slamming event equal to the wave encounter frequency (see Figure 33). The peak values of the slamming pressure are found at the bow flare region (P4), and then the pressure decreases from the bow front to the backward. As shown in Figure 33, the impact loads in the focused wave are generally larger than those in the regular waves among these gauges, with an increase up to 9% of slamming pressure at the bow flare.
Visual observation of slamming events obtained from the CFD analysis at a wave cycle (t = 17 s-18 s) of the ship in regular and focused wave conditions are shown in Figure 34. The pronounced slamming and severe wave overtopping phenomenon are captured well in the focused wave due to the severe ship responses, whereas they do not appear under regular wave conditions.

Discussion
This paper adopts an efficient two-way FSI model coupled with an OpenFOAM solver and MBDyn to investigate the ship hydroelasticity of an S175 containership in regular and focused wave conditions with forward speeds. A series of validation and verification studies were presented to evidence that the present CFD-DMB method can accurately measure the ship responses and the peak VBM loads in waves, which is of prime importance to access ship hydrodynamic and hydroelasticity problems. This well-validated FSI framework is capable of extending its applications for various types of ships at any length. However, a dedicated beam frame needs to be designed for multihull marine structures, e.g., catamaran and trimaran, in reference to [60]. The present method can be further extended to study asymmetrical ship motions, e.g., coupled torsional and bending behaviour of a ship operating in oblique waves. The limitations should be illustrated here: the Euler-Bernoulli beam is not applicable for modelling large open-section bulk carriers or similar cases with ships experiencing huge torsional effects. To do so would require a full three-dimension structural model or a high-order beam, e.g., Vlassov beam.
Determining extreme values of ship responses in irregular waves through FSI simulation based on the Navier-Stokes equations is almost impossible. Here, the focused wave theory is used to represent expected maximum responses from the JONSWAP sea spectrum, i.e., wave events tailored to generate this response. A critical attribute presented by this paper is the comparison of the numerical results, including the ship global motions, ultimate structural strength, slamming and green water loads between focused and regular wave conditions. These data give hints of the significant influence of focused waves on the ship structures, which found that the heave and pitch amplitudes rise by about 15-20% compared to those in regular waves. The VBMs are extracted from the beam profile in the structure solver and compared between the focused and regular wave with time-history data in Section 5.2. These results show that the focused waves have a maximum increase of 8.5% and 9.8% of total hogging and sagging VBMs, respectively, compared to regular waves. The high-frequency components are contained in total VBM signals, which are attributed to the nonlinearity of the response and wave process. In addition, the significant nonlinearities of the response are found as hogging/sagging asymmetry and vibrations from slamming impacts. The slamming impacts between focused waves and regular waves are compared in Section 5.3. We found that the impact pressure signals in focused waves generally show larger and sharper shapes than in regular waves. Moreover, the physical observations of the ship sailing in focused waves clearly show that severe bow slamming and green water on deck occurred; however, it is not noticeable in the regular wave condition.

Conclusions
In this paper, a fully coupled CFD-DMB numerical model has been first developed to study the seakeeping and hydroelastic behaviour of a flexible S175 containership in regular and focused wave conditions. The open-source CFD toolbox OpenFOAM was utilised to model the flow field, and waves were generated by Waves2Foam. A Euler-Bernoulli (EB) beam model, equivalently representing the hull characteristics, was employed in a multibody structure solver, MBDyn, to calculate the structural deformation under excitations. A two-way coupling algorithm was implemented for exchanging data between the fluid and structure solver. The tool is particularly suitable for analysing the hydroelastic effects of slender ships or elastic marine structures in waves containing multiple mechanically connected components. The main conclusions obtained are as follows: 1.
Taking the S175 containership as an example, we have demonstrated the accuracy of this integrated numerical modelling tool in the prediction of dynamic ship motions and slamming wave loads with the hydroelastic effects within a range of wave frequencies. A comprehensive comparison between the present method against the CFD-FEA method and the experimental results show that our CFD-DMB method accurately measures the peak VBM loads, which is persistent in the problems associated with ship hydroelasticity. 2.
The focused wave group was generated based on the modified NewWave theory [31], considering the current effects. The simulated wave elevations were found in favourable agreement with the theory. 3.
The extreme ship motions and hydroelastic responses of the S175 containership were further studied in a focused wave condition to investigate the extreme sea conditions that a ship may experience in a real sea state. Numerical results demonstrated that the flexible ship would experience larger ship motions, vertical bending moments and slamming loads in focused waves than in regular waves, in which case, the design of the hull girder should leave enough of a margin to resist the high loads expected in real sea operations. Our future work will present an experimental study on the flexible ship in focused wave conditions.
Recommendations for future research are briefly outlined below.

1.
The asymmetric loads on the hull girder while the ship is sailing in oblique wave and cross-wave conditions, e.g., coupled horizontal bending moments and torsional moments, could be investigated.

2.
Shallow water wave condition is another research area to take for the investigation of hull girder vibrations induced by the nonlinearity of shallow wave components; however, the ship operational speed will not be as fast as in deep water due to squat effects.

3.
The high-frequency vibration induced by the propeller and machinery is one of the primary sources of causing ship vibrations. A simplified propeller-shaft-ship system should be designed in MBDyn to investigate propeller induced high-frequency vibrations. Data Availability Statement: All data, models, or code generated or used during the study are available from the corresponding author by request.