BEM-FEM Method for the Prediction of Symmetric Hydroelastic Responses of Ships with Forward Speed

,


Introduction
In recent years, the economies of scale lead to increase in the numbers and sizes of large ocean-going vessels.For example, Ultra Large Container Ships (ULCSs) with capacity of more than 20,000 TEU and length and beam exceeding 400 m and 60 m, respectively, is today the norm that challenges the wave load margins introduced by classification rules [1].Because of the slenderness and open deck configuration, these vessels are prone to springing and whipping-induced loads in stochastic seaways [2].Springing is a phenomenon in which wave frequency or its harmonics are able to resonate at the structural natural frequency and may lead to fatigue failure.Whipping results from slamming of ships which causes transient dynamic loading along the hull girder.It is important for structural design as it imparts impact loads leading to fatigue failures [3,4].
Hydroelastic idealisations used for the prediction of springing and whipping loads can be carried out in the frequency or time domains."Frequency-domain" methods gained popularity since late 70s when the basic hydroelasticity theory combining strip theory and Timoshenko beam dynamics was introduced by Bishop and Price [5].e three-dimensional version of this method was introduced by Bishop et al. [6].eir approach combined FEA (for dry analysis) and frequency-domain Green function methods (for the wet analysis).To date, these methods have been applied successfully to a variety of merchant and naval ships, and their capability to simulate symmetric, antisymmetric, and asymmetric dynamic behaviour of ships in waves and for use in design has been widely demonstrated in the literature (e.g., [7][8][9][10][11][12][13][14]).
e concept of frequency-domain modal analysis has been successfully applied not only to describe steady-state responses but also to describe the behaviour of symmetric transient responses due to slamming in irregular seaways [15][16][17].Yet, comparisons between 2D time and frequency-domain techniques demonstrated that the effects of nonlinearities become particularly important at higher speeds and for ships with large flare [18].Since the application of Reynolds-averaged Navier-Stokes computational fluid dynamics (RANS CFD) hydroelastic methods remains computationally time consuming (e.g., [19,20], potential flow time-domain hydroelastic methods remain preferable for the prediction of the influence of whipping and springing loads for ship design and assessment. "Time-domain" potential flow hydroelastic methods can be divided in two categories, namely, (a) the expanded Cummins' equation method and (b) direct simulation methods.In direct simulation methods, structural dynamics are idealised by finite elements and flexible fluid structure interaction (FFSI) is enabled by a boundary element method (BEM) on the interfacing boundaries (e.g., [21][22][23][24][25][26][27][28][29][30][31][32]).Although the direct coupling can be beneficial for strongly nonlinear problems, its applications are computationally expensive.On the other hand, the expansion of Cummins' equation method makes use of "impulse response functions (IRFs)," Fourier transformation, and structural dynamics and is therefore considered simpler in terms of computational modelling and efficiency (e.g., [33][34][35][36][37][38][39]).
Recently, Pal et al. [40,41] used the 3D time-domain panel method based on time-domain free surface Green's function.e ship structure was modelled as an Euler beam.
is work demonstrated that direct coupling hydroelastic methods can be more beneficial for strongly nonlinear problems as they allow for implicit implementation of wet modes and the easier inclusion of hydrodynamic nonlinearities.Weakly nonlinear hydroelasticity methods have also been presented by  and Jiao et al. [39].e former used a Rankine panel method to idealise nonlinear hydrodynamics and a simplified finite element for Vlasov beam dynamics.e latter accounted for the influence of Froude-Krylov forces while radiation-diffraction forces were considered linear.Along these lines, it is now understood that Rankine panel or free surface time-domain panel methods are more suitable for inclusion of nonlinear radiation-diffraction forces.Segmented model tests are widely adopted for the measurement of ship hydroelastic motions and loads (e.g., [33,[42][43][44][45][46]). Notwithstanding this, results from model tests remain uncertain and generally unambiguous [47].
Building up from the work by Pal et al. [40], the originality and focus of this work is on exploring the adequacy of time-domain physical assumptions to model forward speed effects and exploring the influence of Timoshenko beam dynamics (rotary inertia and shear deformation effects) on predicting symmetric flexible dynamic response in waves.As explained in Section 2 of the paper, the formulation of the hydrodynamic problem is based on the model of Datta and Sen [48]; forward ship speed effects are idealised by a space-state function and modal actions are accounted for by Timoshenko beam structural dynamics.Flexible fluid structure interaction (FFSI) coupling is enabled by a body boundary condition, and a direct integration Newmark-β scheme is used to obtain symmetric dynamic responses.Results from a time-domain hydroelastic method that are compared for three different modern container ship hull forms and the experiments of Rajendran and Guedes Soares [33] are presented in Section 3. e influence of structural flexibility on ship responses is discussed in Section 4 and conclusions are drawn in Section 5.

Theory
Figure 1 presents the case of a flexible ship progressing with uniform forward speed U. e origin Oxyz of the vessel lies in the Earth fixed coordinate system in way of the centre of gravity G(X s , Y s , Z s ).
e mathematical idealisation presented in the following sections comprises two parts, namely, "structural dynamic" and "hydrodynamic" idealisations.In the former, FEA-based Timoshenko beam structural dynamics are used to idealise the structural part of the FFSI problem.e boundary integral equation is then formulated using time-domain free surface Green's function, and FFSI coupling is enabled by a direct integration Newmark-β scheme [49].

Structural Dynamic Idealisation.
e structure is assumed to behave as a slender nonuniform Timoshenko beam accounting for shear deformation and rotary inertia effects (see Figure 2) [50].Transverse shear is assumed to be constant over the cross section.e rotation about the y axis is denoted by an independent function, namely, ψ.
e dynamical behaviour of the ship hull due to hydrodynamic external force can be described by governing differential equations.In pure bending, the cross section maintains orthogonality, but in this work, the net slope of the natural axis is presented in terms of both flexure and the shear strain.In this formulation, w and ψ represent independent field variables idealising the transverse deflection of natural fibre and the angle of flexure, respectively.erefore, the shear strains (ε) and stresses (σ) are, respectively, defined as ε � −z and σ � −Ez.e vertical bending moment (M) and shear force (v) are defined as M � EI and v � GA c k s (ψ + dw/dx).In these formulations, G presents the modulus of rigidity (shear modulus) and A s is the effective shear area.As the shear stress is not uniformly distributed over the cross section of the beam, considering uniform area overestimates the stresses.In an actual flexing scenario, the shear force will be lower and the effective shear area is defined as A s � A c k s where A c is the cross-sectional area and k s is the shear correction factor with values generally taken as 1.2 (i.e., approximately 80% of the area if considered effective in shear) for classic ship like beam idealisations [5,51].
e ordinary differential equations expressing the shear forces and bending moments are defined, respectively, by ( 1) and (2) as 2 Shock and Vibration In equation (1), f represents the distributed load over the sectional element of the ship modelled as a beam.e weak form of (1) and (2) can be obtained by introducing two weight functions, namely, transverse deflection w and rotational function ψ, as independent unknowns.A "Gaussian quadrature method (G-Q)" can be applied to evaluate the weak form equations as follows: e weak form of the structural dynamic equations (3) and ( 4) is obtained by introducing two weight function α b and β S ; h e is the length of 1D finite element where x A an d x B are the two consecutive nodes of a beam-like element; Q e 1 and Q e 3 represent elemental shear forces; Q e 2 and Q e 4 denote elemental bending moments.To avoid the shear locking [52], the equations are equally interpolated by a "reduced integration element (RIE)."Accordingly, the elemental static equation can be written as In ( 5), the global stiffness (K ij ) operators are defined as  Shock and Vibration 3 x A GA c k s dψ (1)   i dx dψ (1) e shear stiffness is computed with one Gauss point.e mass matrix is defined as the mass matrices for shear and bending criteria, respectively.erefore, to determine the response history of the beam element, the global finite element equation for dry analysis is expressed as In this equation of motion, the responses of the structure are demonstrated in the form of displacement (x), velocity ( _ x), and acceleration ( € x); [M] and [K] represent the global mass and stiffness matrices, [C] is the damping matrix, and F ext   is the global external force vector.e formulation presented assumes Rayleigh damping [53] where e value of stiffness proportional coefficient (σ d ) and mass proportional coefficient (ς d ) is taken as presented by Liu et al. [54].

Hydrodynamic Idealisation.
In (7), the external force F ext   vector represents the wave induced loads and moments which may be calculated by solving the hydrodynamic problem.In this study, a 3D time-domain panel method based on transient free surface Green's function has been employed.e theory is formulated based on Earth fixed coordinate system.e complete formulation of the hydrodynamic problem for rigid body case is presented in multiple sources (e.g., [48,55]).erefore, the background to the method is only briefly presented here with the aim to explain its relevance for the solution of the hydroelasticity problem.In Figure 3, Ω represents a fluid domain of interest which is surrounded by the free surface S F , bottom surface S b , wetted body surface S 0 , and the surfaces S −∞ , S ∞ at infinity.e total velocity potential ϕ T (p, t) is defined as where ϕ I (p, t) is the velocity potential of incident waves at any time instant (t) at any arbitrary point p on the fluid domain Ω and ϕ(p, t) is the disturbed potential that consists of radiation and diffraction components.Within the context of linear potential theory, solution for ϕ(p, t) can be obtained by solving following governing equation: Using linear free surface boundary condition, body boundary condition, bottom boundary condition, and radiation condition as follows: within the context of rigid body ship dynamics, we can assume that the normal velocity of the body is taken as a function of time.However, if we account for the influence of hydroelasticity, the velocity of the body should be idealised by a state-space function.Considering that the formulation is based on Earth fixed coordinate system, consideration of the so called "m terms" [26] is not required.To solve the hydrodynamic problem, transient free surface Green's function is used.e expression for ϕ(p, t) becomes 4 Shock and Vibration e detailed expression for the derivation of ( 11) is given in Lin and Yue [55]; G o presents the nonlinear part (Rankine part) of Green's function and G f t is the linear part (regular part) of Green's function [56].
is solution produces scattering of velocity potential ϕ(p, t) over the surface of the hull; v N and v n are normal to the water line Γ and the wetted surface, respectively; σ(p, t) represents the strength of the source G; t and τ represent the present and past time increments used in the convolution time integral.e time dependent contribution comes from the second and third integrals that calculate the effect that comes from the disturbance caused in previous time steps.e third integral arises due to the forward speed effects with terms "t" and "τ" denoting the present and past instants in time.
e computation of the velocity potential term ϕ(p, t) leads to prediction of dynamic pressure P(p, t) defined as In ( 12), the density of the fluid is termed as ρ 0 .For zero speed, the contribution of quadratic term can be ignored.However, when forward speed effects are accounted for, consideration of quadratic terms is essential.Accordingly, the method is not fully linear and fundamentally different from other time-domain formulations (e.g., [26,39]).Solving (12) gives the distribution of the pressure at any arbitrary point at any instant of time t.If we assume that P(X, t) represents the collection of total dynamic pressure at any point on the sectional curve C S x for an arbitrary vertical section S x , the sectional force F S x is defined as where (n) denotes the generalized normal on a section increment (dS X ). e restoring force F restoring i (t) at an i th panel of a surface increment ds i subject to normal n i may expressed as where Z i is the water head and δ denotes the leading order variation.Lower-order panel methods assume that variation over a hydrodynamic panel remains constant.As δ varies over the surface, this leading order variation might not be very effective, and therefore Kim et al. [26] proposed the use of a higher-order Rankine panel method.Notwithstanding this, the stability of the numerical idealisation can also be ensured by reducing the size of panels (i.e., using optimum number of fluid domain discretisation) when using a classic Green function method (e.g., [57]).

Flexible Fluid Structure Interaction (FFSI)
. Equation ( 7) may be solved by Newmark [49] time integration method.Accordingly, the nodal displacement vector x { } at any time instant t is obtained from the equation e stability and accuracy of the present numerical scheme during a time increment Δt is dependent on the parameters α and δ. ese two parameters also describe the variation of the acceleration over a time step.e numerical values of α and δ in the above equation are in general taken as ¼ and ½, respectively (Kim et al., 2013).e time variations of the acceleration and velocity derivatives are defined as

Shock and Vibration
Accordingly, the velocity and displacement at each section may be obtained from equations ( 15)- (17).By back substitution in equation ( 14), the boundary condition for the solution of the fluid problem is obtained and FFSI is enabled.

Results
Table 1 outlines the principal particulars of three container ship hulls that have been assessed using the method presented in this publication.e hydrodynamic idealisations for each design and the body plans are shown in Figure 4.
e weight and flexural rigidity distributions for all containerships are depicted in Figure 5. Ship 1 is the well-known S175 container ship [46].Ship 2 is a modern large container ship (LC) design studied under previous ISSC-ITTC benchmarks [3,33,58].Ship 3 is a modern ultra-large container ship (ULC) [59].
e results presented in this paper are given in the nondimensional forms presented in Table 2. Numerical results based on the theory presented in Section 2 are named as TDFlex-T.e computations based on Datta and Soares [41] model are denoted as TDFlex-E, and the results from time-domain rigid body approach of Sengupta et al. [60] are referenced as TDRigid.

Validation Study. Figures 5(a)-5(d) demonstrate the time history of nondimensional vertical bending moments
(VBMs) and shear forces (VSFs) of S175 and ULC at the wave matching region (i.e., λ/L � 1) for forward speed F n � 0.25.e results presented converge to steady state.
is confirms the numerical stability of the solution.In Tables 3-5, dry natural frequencies for the different modes of S175, LC, and ULC are compared against the Euler-Bernoulli models of Rajendran et al. [59] and Rajendran and Guedes Soares [33] as well as the superposition method of Wu and Hermundstad [18].Variations in the natural frequencies can be attributed to the Timoshenko beam idealisation introduced in Section 2.
In Figure 6, the VBM transfer function for S175 hull at F n � 0.25 is plotted and compared with the numerical method proposed by Datta and Guedes Soares [41] that has been based on an Euler beam model.the experimental results published by Iijima et al. [61].e differences in the magnitude of the VBM value at the wave matching region (λ/L � 1) could be attributed to limitations of the linear hydrodynamic theory assumptions and structural idealisations.Notwithstanding this, it is thought that the approach presented captures very well the overall trend of the dynamic response.
Comparisons of results for amidships symmetric responses against [33] are depicted in Figure 7. Figure 7(a) presents a comparison of nondimensional displacement against nondimensional frequency.As shown in Figure 7(b), the magnitude of the VBM matches well with the other published results.Comparison of VBM along the hull for LC is demonstrated in Figure 8 at F n � 0.12.Once again, the results compare well with experiments and the coupled BEM-FEM model given by Datta and Guedes Soares [41].
e VBM at amidships the ULC for F n = 0.135 is shown in Figure 9. Results are compared with experiments by Rajendran et al. [59] and TDFlex-E.Both results underestimate the maximum VBM.However, the overall trend agrees very well with experimental values.It may be therefore concluded that the present numerical model is stable and consistent and captures the symmetric responses with reasonable accuracy.

Parametric Study
In this section, a parametric study is carried out to observe the effect of design parameters in symmetric hydroelastic responses such as vertical displacement, VSF, and VBM. e parametric study presented utilises S175 and ULC hulls.Special focus is attributed to the influence of structural rigidity, forward speed, ship length, and slenderness on dynamic response.

e Influence of Structural Rigidity.
e vertical displacement RAO at midships and transfer functions of VSF at + L/4position from the aft of the ULC were compared between the flexible and rigid structures for F n � 0.135.Figure 10 shows the comparison of the displacement response amplitude operator (RAO) for rigid and flexible structures over the nondimensional frequency ω ��� L/g  .Observed differences demonstrate that rigid body assumptions overestimate the structural responses.As shown in Figure 11, near the resonating frequency, the ship is likely to experience larger shear forces, which can only be captured if structural flexibility is taken into account.While TDFlex-T is capable of capturing this phenomenon, TDRigid fails to capture the same phenomenon.is implies the importance of capturing the influence of structural flexibility for slender vessels.It is observed that the peak symmetric loads drop rapidly with respect to the frequency.is implies that the chances of damage of the structure are higher at lower Froude number at lower frequency zone.

6
Shock and Vibration     14). is could be attributed to the influence of resonance phenomena.On the other hand, the VSF values for the S175 hull decrease slowly and do not show any secondary peak in the response (see Figures 14-16). is nature of the transfer function of the VSF values for the S175 hull is similar to the rigid body case.It may therefore be concluded that hydroelasticity effects may be more prominent for more slender vessels.

e Influence of Structural Flexibility on Hydrodynamic Pressure Distributions.
e influence of flexibility on local pressure variation was evaluated at the free surface and the hull bottom in way of the bow, stern, and amidships for F n � 0.25 and λ/L � 1 (see Figure 17). is is because these conditions were shown to be critical in terms of flexible fluid structure interactions (see Section 4).
e variations of nondimensional hydrodynamic pressures over time are shown in Figures 18-20.Figures 18(a          Shock and Vibration

Conclusions
is paper presented a direct hydroelastic analysis method that combines a time-domain Green function with Timoshenko beam structural dynamics for the prediction of symmetric ship responses in waves.
e hydrodynamic model is based on the linear 3D time-domain panel method where the Timoshenko beam model is adopted for the structural solution.Results were validated by comparisons against a range of available experimental data.A parametric study for three different containership hull forms demonstrated that hull slenderness and forward speed may influence symmetric flexible ship responses.Future research will focus on developing a fully nonlinear hydroelastic method for the prediction of extreme sea loads (e.g., slamming, green water on decks, and so on) on hull forms of low rigidity progressing at a medium to high speed in irregular waves.

A:
Wave amplitude (m) E: Young's modulus (N/m 2 ) I: Moment of inertia (kg-m 2 ) F n : Froude number f: Sectional load (N/m) G: Shear modulus (N/m 2 ) G(x s y s z s ): Ship centre of gravity g: Acceleration of gravity (m/sec 2 ) K ij : Global stiffness factor k s : Shear correction factor L: Hull length (m) M: Vertical bending moment (N-m) Oxyz: Vessel Earth fixed coordinate system in way of the centre of gravity p: Field point q: Source point S 0 : Mean wetted surface T: Time period (sec) t: Time (sec) U: Forward speed (m/s) V: Vertical shear force (N) w: Transverse deflection (m) x 3 : Vertical displacement (m) Greek symbols α b , β s : Residual weight functions λ: Wavelength (m) ρ: Structural density (kg/m 3 ) ρ o : Water density (kg/m 3 ) ψ: Bending rotation (rad) ω: Wave frequency (Hz).

Figure 2 :
Figure 2: Section of a beam representing shear forces (V) and moments (M) at edges.

▽ 2 Figure 3 :
Figure 3: Hydrodynamic boundary domain and associated conditions (S denotes surface of relevance to boundary conditions).

14
Shock and Vibration amidships (see Figure17; location points 2 and 5) and astern (see Figure17; location points 3 and 6).e variation of the pressure along the ship hull for a particular time instant is depicted in Figures21(a) and 21(b).e same plots display that the pressure variation is more visible for the case of ULC.It may be therefore concluded that as the ship length increases, the influence of hull flexibility on hydrodynamic pressures becomes prominent.

Table 1 :
General particulars of different container ship models.

Table 3 :
Dry natural frequencies of the S175 for vertical flexural vibrations.

Table 4 :
Dry natural frequencies of the LCS for vertical flexural vibrations.

Table 5 :
Dry natural frequencies of the ULCS for vertical flexural vibrations.