# A Fully Coupled Time-Domain BEM-FEM Method for the Prediction of Symmetric Hydroelastic Responses of Ships with Forward Speed

Authors: Tushar Kanti Show; Spyros Hirdaris; Ranadev Datta
Journal: Shock and Vibration
 (2022)

---

## Abstract

This paper presents a direct time-domain method for the prediction of symmetric hydroelastic responses of ships progressing with forward speed in small amplitude waves. A transient time-domain free surface Green function is used for the idealisation of the seakeeping problem using an Earth fixed coordinate system. Free surface ship hydrodynamics are idealised in the time domain by a Green function, and forward speed effects are idealised by a space-state model. 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. The method is validated against available published numerical and experimental results. A parametric study for different container ship hull forms confirms that (i) forward speed effects should be taken under consideration as far as practically possible and (ii) hull flexibility effects accounting for hull shear deformation and rotary inertia are more notable for slender hull forms.


---

## Body

## 1. 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]. The three-dimensional version of this method was introduced by Bishop et al. [6]. Their 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–14]).The 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–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–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–39]).Recently, Pal et al. [40, 41] used the 3D time-domain panel method based on time-domain free surface Green’s function. The ship structure was modelled as an Euler beam. This 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 Kim et al. [26–28] and Jiao et al. [39]. The former used a Rankine panel method to idealise nonlinear hydrodynamics and a simplified finite element for Vlasov beam dynamics. The 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–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. The influence of structural flexibility on ship responses is discussed in Section 4 and conclusions are drawn in Section 5.


## 2. Theory

Figure1 presents the case of a flexible ship progressing with uniform forward speed U. The origin Oxyz of the vessel lies in the Earth fixed coordinate system in way of the centre of gravity GXs,Ys,Zs. The 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. The 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].Figure 1
Schematic diagram of coordinate systems.


### 2.1. Structural Dynamic Idealisation

The structure is assumed to behave as a slender nonuniform Timoshenko beam accounting for shear deformation and rotary inertia effects (see Figure2) [50]. Transverse shear is assumed to be constant over the cross section. The rotation about the y axis is denoted by an independent function, namely, ψ. The 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. Therefore, the shear strains (ε) and stresses (σ) are, respectively, defined as ε=−z and σ=−Ez. The vertical bending moment (M) and shear force v are defined as M=EI and v=GAcksψ+dw/dx. In these formulations, G presents the modulus of rigidity (shear modulus) and As 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 As=Acks where Ac is the cross-sectional area and ks 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].Figure 2
Section of a beam representing shear forces (V) and moments (M) at edges.
(a)(b)(c)The ordinary differential equations expressing the shear forces and bending moments are defined, respectively, by (1) and (2) as(1)ddxGAcksψ+dwdx+f=0,(2)ddxEIdψdx−GAcksψ+dwdx=0.In equation (1), f represents the distributed load over the sectional element of the ship modelled as a beam. The 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:(3)∫0hedwdxdαbdxGAcksdx+ψdαbdxGAcksdx=∫0hefαbdx+αbheQ1e+Q3eαbhe,(4)∫0heβsdwdxGAcksdx+∫0hedβsdxdψdxEI+βsψGAcksdx=βsxAQ2e+Q4eβshe.The weak form of the structural dynamic equations (3) and (4) is obtained by introducing two weight function αb and βS; he is the length of 1D finite element where xAandxB are the two consecutive nodes of a beam-like element; Q1e and Q3e represent elemental shear forces; Q2e and Q4e 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(5)K11K12K21K22wψ=f1f2.In (5), the global stiffness (Kij) operators are defined as(6)Kij11=∫xAxBGAcksdψi1dxdψj1dxdx,i=1,2,3,……nj=1,2,3,……n,Kij12=∫xAxBGAcksdψi1dxψj2dx,i=1,2,3,……nj=1,2,3,……n,Kij21=∫xAxBGAcksdψi1dxψj2dx,i=1,2,3,……nj=1,2,3,……n,Kij22=∫xAxBGAcksψi2ψj2dx+∫xAxBEIdψi2dxdψj2dxdx,i=1,2,3,……nj=1,2,3,……n.The shear stiffness is computed with one Gauss point. The mass matrix is defined asmsb=ms+mb, where ms and mb represent the mass matrices for shear and bending criteria, respectively. Therefore, to determine the response history of the beam element, the global finite element equation for dry analysis is expressed as(7)Mx˙+Cx˙+Kx=Fext.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 Fext is the global external force vector. The formulation presented assumes Rayleigh damping [53] where C=σdK+ςdM. The value of stiffness proportional coefficient (σd) and mass proportional coefficient (ςd) is taken as presented by Liu et al. [54].


### 2.2. Hydrodynamic Idealisation

In (7), the external force Fext 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. The theory is formulated based on Earth fixed coordinate system. The complete formulation of the hydrodynamic problem for rigid body case is presented in multiple sources (e.g., [48, 55]). Therefore, 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 SF, bottom surface Sb, wetted body surface S0, and the surfaces S−∞, S∞ at infinity.Figure 3
Hydrodynamic boundary domain and associated conditions (S denotes surface of relevance to boundary conditions).The total velocity potentialϕTp,t is defined as(8)ϕTp,t=ϕIp,t+ϕp,t,where ϕIp,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:(9)∇2ϕp,t=0,forp∈Ω.Using linear free surface boundary condition, body boundary condition, bottom boundary condition, and radiation condition as follows:(10)∂2ϕ∂t2+g∂ϕ∂z=0atz=0,∂ϕ∂n=vnx,t−∂ϕI∂nonS0,∂ϕ∂z=0onSbϕ,∂ϕ∂t=0atS∞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. The expression for ϕp,t becomes(11)ϕp,t=−14π∬S0σq,tGop,qdS+∫0tdτ∬S0σq,tGtfp,q;t−τdS−∫0τdτ1g∫Γσq,tGtfp,q;t−τvNvndL.The detailed expression for the derivation of (11) is given in Lin and Yue [55]; Go presents the nonlinear part (Rankine part) of Green’s function and Gtf is the linear part (regular part) of Green’s function [56]. This solution produces scattering of velocity potential ϕp,t over the surface of the hull; vN and vn 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. The time dependent contribution comes from the second and third integrals that calculate the effect that comes from the disturbance caused in previous time steps. The third integral arises due to the forward speed effects with terms “t” and “τ” denoting the present and past instants in time.The computation of the velocity potential termϕp,t leads to prediction of dynamic pressure Pp,t defined as(12)Pp,t=−ρ0ϕtp,t−Uϕxp,t+ϕtIp,t+12∇ϕp,t+∇ϕIp,t2.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 PX¯,t represents the collection of total dynamic pressure at any point on the sectional curve CSx for an arbitrary vertical section Sx, the sectional force FSx is defined as(13)FSx=∫CSxPX¯,t.n⌢.dCSx,where n¯ denotes the generalized normal on a section increment (dSX).The restoring forceFirestoringt at an ith panel of a surface increment dsi subject to normal n¯i may expressed as(14)Firestoringt=−ρ0gZ¯iδn¯i.dsi−ρ0gδZ¯in¯i.dsi,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]).


### 2.3. 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(15)MαΔt2+δαΔtC+Kxt=Ft+M1αΔt2xt−Δt+1αΔtx˙t−Δt+12α−1x¨t−Δt+CδαΔtxt−Δt+δα−1x˙t−Δt+Δt2δα−2x¨˙t−Δt.The stability and accuracy of the present numerical scheme during a time incrementΔt is dependent on the parameters α and δ. These two parameters also describe the variation of the acceleration over a time step. The numerical values of α and δ in the above equation are in general taken as ¼ and ½, respectively (Kim et al., 2013). The time variations of the acceleration and velocity derivatives are defined as(16)x¨t=1αΔt2xt−xt−Δt−1αΔtx˙t−Δt−12α−1x¨t−Δt,(17)x˙t=x˙t−Δt+Δt1−δx¨t−Δt+δx¨t.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.


## 2.1. Structural Dynamic Idealisation

The structure is assumed to behave as a slender nonuniform Timoshenko beam accounting for shear deformation and rotary inertia effects (see Figure2) [50]. Transverse shear is assumed to be constant over the cross section. The rotation about the y axis is denoted by an independent function, namely, ψ. The 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. Therefore, the shear strains (ε) and stresses (σ) are, respectively, defined as ε=−z and σ=−Ez. The vertical bending moment (M) and shear force v are defined as M=EI and v=GAcksψ+dw/dx. In these formulations, G presents the modulus of rigidity (shear modulus) and As 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 As=Acks where Ac is the cross-sectional area and ks 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].Figure 2
Section of a beam representing shear forces (V) and moments (M) at edges.
(a)(b)(c)The ordinary differential equations expressing the shear forces and bending moments are defined, respectively, by (1) and (2) as(1)ddxGAcksψ+dwdx+f=0,(2)ddxEIdψdx−GAcksψ+dwdx=0.In equation (1), f represents the distributed load over the sectional element of the ship modelled as a beam. The 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:(3)∫0hedwdxdαbdxGAcksdx+ψdαbdxGAcksdx=∫0hefαbdx+αbheQ1e+Q3eαbhe,(4)∫0heβsdwdxGAcksdx+∫0hedβsdxdψdxEI+βsψGAcksdx=βsxAQ2e+Q4eβshe.The weak form of the structural dynamic equations (3) and (4) is obtained by introducing two weight function αb and βS; he is the length of 1D finite element where xAandxB are the two consecutive nodes of a beam-like element; Q1e and Q3e represent elemental shear forces; Q2e and Q4e 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(5)K11K12K21K22wψ=f1f2.In (5), the global stiffness (Kij) operators are defined as(6)Kij11=∫xAxBGAcksdψi1dxdψj1dxdx,i=1,2,3,……nj=1,2,3,……n,Kij12=∫xAxBGAcksdψi1dxψj2dx,i=1,2,3,……nj=1,2,3,……n,Kij21=∫xAxBGAcksdψi1dxψj2dx,i=1,2,3,……nj=1,2,3,……n,Kij22=∫xAxBGAcksψi2ψj2dx+∫xAxBEIdψi2dxdψj2dxdx,i=1,2,3,……nj=1,2,3,……n.The shear stiffness is computed with one Gauss point. The mass matrix is defined asmsb=ms+mb, where ms and mb represent the mass matrices for shear and bending criteria, respectively. Therefore, to determine the response history of the beam element, the global finite element equation for dry analysis is expressed as(7)Mx˙+Cx˙+Kx=Fext.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 Fext is the global external force vector. The formulation presented assumes Rayleigh damping [53] where C=σdK+ςdM. The value of stiffness proportional coefficient (σd) and mass proportional coefficient (ςd) is taken as presented by Liu et al. [54].


## 2.2. Hydrodynamic Idealisation

In (7), the external force Fext 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. The theory is formulated based on Earth fixed coordinate system. The complete formulation of the hydrodynamic problem for rigid body case is presented in multiple sources (e.g., [48, 55]). Therefore, 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 SF, bottom surface Sb, wetted body surface S0, and the surfaces S−∞, S∞ at infinity.Figure 3
Hydrodynamic boundary domain and associated conditions (S denotes surface of relevance to boundary conditions).The total velocity potentialϕTp,t is defined as(8)ϕTp,t=ϕIp,t+ϕp,t,where ϕIp,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:(9)∇2ϕp,t=0,forp∈Ω.Using linear free surface boundary condition, body boundary condition, bottom boundary condition, and radiation condition as follows:(10)∂2ϕ∂t2+g∂ϕ∂z=0atz=0,∂ϕ∂n=vnx,t−∂ϕI∂nonS0,∂ϕ∂z=0onSbϕ,∂ϕ∂t=0atS∞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. The expression for ϕp,t becomes(11)ϕp,t=−14π∬S0σq,tGop,qdS+∫0tdτ∬S0σq,tGtfp,q;t−τdS−∫0τdτ1g∫Γσq,tGtfp,q;t−τvNvndL.The detailed expression for the derivation of (11) is given in Lin and Yue [55]; Go presents the nonlinear part (Rankine part) of Green’s function and Gtf is the linear part (regular part) of Green’s function [56]. This solution produces scattering of velocity potential ϕp,t over the surface of the hull; vN and vn 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. The time dependent contribution comes from the second and third integrals that calculate the effect that comes from the disturbance caused in previous time steps. The third integral arises due to the forward speed effects with terms “t” and “τ” denoting the present and past instants in time.The computation of the velocity potential termϕp,t leads to prediction of dynamic pressure Pp,t defined as(12)Pp,t=−ρ0ϕtp,t−Uϕxp,t+ϕtIp,t+12∇ϕp,t+∇ϕIp,t2.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 PX¯,t represents the collection of total dynamic pressure at any point on the sectional curve CSx for an arbitrary vertical section Sx, the sectional force FSx is defined as(13)FSx=∫CSxPX¯,t.n⌢.dCSx,where n¯ denotes the generalized normal on a section increment (dSX).The restoring forceFirestoringt at an ith panel of a surface increment dsi subject to normal n¯i may expressed as(14)Firestoringt=−ρ0gZ¯iδn¯i.dsi−ρ0gδZ¯in¯i.dsi,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]).


## 2.3. 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(15)MαΔt2+δαΔtC+Kxt=Ft+M1αΔt2xt−Δt+1αΔtx˙t−Δt+12α−1x¨t−Δt+CδαΔtxt−Δt+δα−1x˙t−Δt+Δt2δα−2x¨˙t−Δt.The stability and accuracy of the present numerical scheme during a time incrementΔt is dependent on the parameters α and δ. These two parameters also describe the variation of the acceleration over a time step. The numerical values of α and δ in the above equation are in general taken as ¼ and ½, respectively (Kim et al., 2013). The time variations of the acceleration and velocity derivatives are defined as(16)x¨t=1αΔt2xt−xt−Δt−1αΔtx˙t−Δt−12α−1x¨t−Δt,(17)x˙t=x˙t−Δt+Δt1−δx¨t−Δt+δx¨t.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.


## 3. Results

Table1 outlines the principal particulars of three container ship hulls that have been assessed using the method presented in this publication. The hydrodynamic idealisations for each design and the body plans are shown in Figure 4. The 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].Table 1
General particulars of different container ship models.
General particularsShip typeNotation (unit)ItemS175LCULCLPP (m)Length between perpendiculars175.0286.6333.43B (m)Beam25.440.042.8D (m)Depth11.024.227.3T (m)Draft9.511.9813.1M (tonne)Total mass24560.085562.71.25e + 05LCG (m)Longitudinal centre of gravity2.784.94.7kxx (m)Radii of gyration along longitudinal (xx); transverse (yy); and vertical (zz) axes10.1614.416.94kyy (m)42.870.14483.9kzz (m)43.011.41284.1Figure 4
Body plan and hydrodynamic idealisations for different container ship models. (a) Ship 1: S175. (b) Ship 2: LC. (c) Ship 3: ULC. (d) Mass distributions. (e) Flexural rigidity distributions.
(a)(b)(c)(d)(e)Figure 5
Symmetric load time histories for ULC and S175 container ships (Fn=0.25, λ/L=1. (a) VBM time history, S175. (b) VSF time history, S175. (c) VBM time history, ULC. (d) VSF time history, ULC.
(a)(b)(c)(d)The results presented in this paper are given in the nondimensional forms presented in Table2. Numerical results based on the theory presented in Section 2 are named as TDFlex-T. The 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.Table 2
Nondimensional dynamic response parameters.
ParametersNondimensional formatFrequency (Hz)ωL/g, λ/LTime (sec)t/TVertical response (m)x3/ALength along the horizontal axisx/LVSF (N/m2)V/ρ0gAL2VBM (N-m)M/ρ0gAL3


### 3.1. Validation Study

Figures5(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 Fn=0.25. The results presented converge to steady state. This 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.Table 3
Dry natural frequencies of the S175 for vertical flexural vibrations.
Mode no.TDFlex-T (Hz)Wu and Hermundstad [18] (Hz)Difference (%)10.30020.3024−0.73720.72990.7669−0.95031.39811.4006−0.170Table 4
Dry natural frequencies of the LCS for vertical flexural vibrations.
Mode no.TDFlex-T (Hz)Rajendran and Guedes Soares [33] (Hz)Difference (%)10.6250.64−2.34321.6031.64−2.25633.0863.11−0.772Table 5
Dry natural frequencies of the ULCS for vertical flexural vibrations.
Mode no.TDFlex-T (Hz)Rajendran et al. [59] (Hz)Difference (%)30.7170.72−0.41641.8641.89−1.3753.5273.64−3.10In Figure6, the VBM transfer function for S175 hull at Fn=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]. The 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.Figure 6
VBM transfer function at midship of S175 (Fn=0.25).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 Fn=0.12. Once again, the results compare well with experiments and the coupled BEM-FEM model given by Datta and Guedes Soares [41].Figure 7
Symmetric responses, LC (atFn = 0). (a) Vertical displacement. (b) VBM.
(a)(b)Figure 8
VBMs along the length for LC (Fn=0.12 and λ/L=1.06).The VBM at amidships the ULC forFn = 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.Figure 9
Comparison of rigid body and hydroelastic RAO of vertical deflection for the ULC at L/4.


## 3.1. Validation Study

Figures5(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 Fn=0.25. The results presented converge to steady state. This 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.Table 3
Dry natural frequencies of the S175 for vertical flexural vibrations.
Mode no.TDFlex-T (Hz)Wu and Hermundstad [18] (Hz)Difference (%)10.30020.3024−0.73720.72990.7669−0.95031.39811.4006−0.170Table 4
Dry natural frequencies of the LCS for vertical flexural vibrations.
Mode no.TDFlex-T (Hz)Rajendran and Guedes Soares [33] (Hz)Difference (%)10.6250.64−2.34321.6031.64−2.25633.0863.11−0.772Table 5
Dry natural frequencies of the ULCS for vertical flexural vibrations.
Mode no.TDFlex-T (Hz)Rajendran et al. [59] (Hz)Difference (%)30.7170.72−0.41641.8641.89−1.3753.5273.64−3.10In Figure6, the VBM transfer function for S175 hull at Fn=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]. The 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.Figure 6
VBM transfer function at midship of S175 (Fn=0.25).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 Fn=0.12. Once again, the results compare well with experiments and the coupled BEM-FEM model given by Datta and Guedes Soares [41].Figure 7
Symmetric responses, LC (atFn = 0). (a) Vertical displacement. (b) VBM.
(a)(b)Figure 8
VBMs along the length for LC (Fn=0.12 and λ/L=1.06).The VBM at amidships the ULC forFn = 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.Figure 9
Comparison of rigid body and hydroelastic RAO of vertical deflection for the ULC at L/4.


## 4. 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. The 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.


### 4.1. The Influence of Structural Rigidity

The 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 Fn=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. This implies the importance of capturing the influence of structural flexibility for slender vessels.Figure 10
RAOs for vertical displacement for flexible body and rigid body, ULC.Figure 11
Comparison of rigid body and hydroelastic RAO of VSF for the ULC at L/4.


### 4.2. The Influence of Forward Speed

Figures12(a), 12(b), 13(a), and 13(b) represent the VBM and VSF for S175 hull and the ULC hull, respectively, for different Froude numbers (Fn=0.135,0.2, and 0.25). The VBM and VSF results are plotted against the nondimensional frequency ωL/g. It is observed that the peak symmetric loads drop rapidly with respect to the frequency. This implies that the chances of damage of the structure are higher at lower Froude number at lower frequency zone.Figure 12
(a) Transfer function for VBM at amidships for different forward speeds of S175. (b) Transfer function for VSF at amidships for different forward speeds of S175.
(a)(b)Figure 13
(a) Transfer function for VBM at amidships for different forward speeds of ULC. (b) Transfer function for VSF at amidships for different forward speeds of ULC.
(a)(b)


### 4.3. The Influence of Ship Length

VSF values for the ULC decrease rapidly at aroundωL/g=2 and then increase (see Figure 14). This 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). This 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.Figure 14
Shear force RAO comparison between S175 and the ULC (Fn = 0.135).Figure 15
Shear force RAO comparison between S175 and the ULC (Fn = 0.2).Figure 16
Shear force RAO comparison between S175 and the ULC (Fn = 0.25).


### 4.4. The Influence of Structural Flexibility on Hydrodynamic Pressure Distributions

The 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 forFn=0.25 and λ/L=1 (see Figure 17). This is because these conditions were shown to be critical in terms of flexible fluid structure interactions (see Section 4).Figure 17
Schematic diagram of selection of panels (squares 1–3 highlight free surface locations; 4–6 highlight ship bottom locations).The variations of nondimensional hydrodynamic pressures over time are shown in Figures18–20. Figures 18(a) and 18(b) display the pressure time histories in way of the bow region of both ships (Figure 17; location points 1 and 4). It is shown that the influence of flexibility is significant only for the ULC container ship model. This trend is confirmed by the comparisons shown in Figures 19(a), 19(b), 20(a), 20(b), respectively, corresponding to responses in way of amidships (see Figure 17; location points 2 and 5) and astern (see Figure 17; location points 3 and 6). The variation of the pressure along the ship hull for a particular time instant is depicted in Figures 21(a) and 21(b). The 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.Figure 18
(a) Pressure time history of S175 hull at bow region,Fn=0.25, λ/L=1. (A) Pressure on free surface (see Figure 17; location point 1). (B) Pressure on bottom surface (see Figure 17; location point 4). (b) Pressure time history of ULC at bow region, Fn=0.25, λ/L=1. (A) Pressure at free surface (Figure 17; location point 1). (B) Pressure at bottom surface (Figure 17; location point 4).
(a)(b)Figure 19
(a) Pressure time history of S175 hull at midship region,Fn=0.25, λ/L=1. (A) Pressure at free surface (Figure 17; location point 2). (B) Pressure at bottom surface (Figure 17; location point 5). (b) Pressure time history of ULC at midship region, Fn=0.25, λ/L=1. (A) Pressure at free surface (Figure 17; location point 2). (B) Pressure at bottom surface (Figure 17; location point 5).
(a)(b)Figure 20
(a) Pressure time history of S175 hull at stern region,Fn=0.25, λ/L=1. (A) Pressure at free surface (Figure 17; location point 3). (B) Pressure at bottom surface (Figure 17, location point 6). (b) Pressure time history of S175 hull at stern region, Fn=0.25, λ/L=1. (A) Pressure at free surface (Figure 17; location point 3). (B) Pressure at bottom surface (Figure 17; location point 6).
(a)(b)Figure 21
(a) Contour plot of pressure profile for a particular instant of time of S175 ship,Fn=0.25, λ/L=1. (b) Contour plot of pressure profile for a particular instant of time of ULC, Fn=0.25, λ/L=1.
(a)(b)


## 4.1. The Influence of Structural Rigidity

The 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 Fn=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. This implies the importance of capturing the influence of structural flexibility for slender vessels.Figure 10
RAOs for vertical displacement for flexible body and rigid body, ULC.Figure 11
Comparison of rigid body and hydroelastic RAO of VSF for the ULC at L/4.


## 4.2. The Influence of Forward Speed

Figures12(a), 12(b), 13(a), and 13(b) represent the VBM and VSF for S175 hull and the ULC hull, respectively, for different Froude numbers (Fn=0.135,0.2, and 0.25). The VBM and VSF results are plotted against the nondimensional frequency ωL/g. It is observed that the peak symmetric loads drop rapidly with respect to the frequency. This implies that the chances of damage of the structure are higher at lower Froude number at lower frequency zone.Figure 12
(a) Transfer function for VBM at amidships for different forward speeds of S175. (b) Transfer function for VSF at amidships for different forward speeds of S175.
(a)(b)Figure 13
(a) Transfer function for VBM at amidships for different forward speeds of ULC. (b) Transfer function for VSF at amidships for different forward speeds of ULC.
(a)(b)


## 4.3. The Influence of Ship Length

VSF values for the ULC decrease rapidly at aroundωL/g=2 and then increase (see Figure 14). This 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). This 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.Figure 14
Shear force RAO comparison between S175 and the ULC (Fn = 0.135).Figure 15
Shear force RAO comparison between S175 and the ULC (Fn = 0.2).Figure 16
Shear force RAO comparison between S175 and the ULC (Fn = 0.25).


## 4.4. The Influence of Structural Flexibility on Hydrodynamic Pressure Distributions

The 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 forFn=0.25 and λ/L=1 (see Figure 17). This is because these conditions were shown to be critical in terms of flexible fluid structure interactions (see Section 4).Figure 17
Schematic diagram of selection of panels (squares 1–3 highlight free surface locations; 4–6 highlight ship bottom locations).The variations of nondimensional hydrodynamic pressures over time are shown in Figures18–20. Figures 18(a) and 18(b) display the pressure time histories in way of the bow region of both ships (Figure 17; location points 1 and 4). It is shown that the influence of flexibility is significant only for the ULC container ship model. This trend is confirmed by the comparisons shown in Figures 19(a), 19(b), 20(a), 20(b), respectively, corresponding to responses in way of amidships (see Figure 17; location points 2 and 5) and astern (see Figure 17; location points 3 and 6). The variation of the pressure along the ship hull for a particular time instant is depicted in Figures 21(a) and 21(b). The 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.Figure 18
(a) Pressure time history of S175 hull at bow region,Fn=0.25, λ/L=1. (A) Pressure on free surface (see Figure 17; location point 1). (B) Pressure on bottom surface (see Figure 17; location point 4). (b) Pressure time history of ULC at bow region, Fn=0.25, λ/L=1. (A) Pressure at free surface (Figure 17; location point 1). (B) Pressure at bottom surface (Figure 17; location point 4).
(a)(b)Figure 19
(a) Pressure time history of S175 hull at midship region,Fn=0.25, λ/L=1. (A) Pressure at free surface (Figure 17; location point 2). (B) Pressure at bottom surface (Figure 17; location point 5). (b) Pressure time history of ULC at midship region, Fn=0.25, λ/L=1. (A) Pressure at free surface (Figure 17; location point 2). (B) Pressure at bottom surface (Figure 17; location point 5).
(a)(b)Figure 20
(a) Pressure time history of S175 hull at stern region,Fn=0.25, λ/L=1. (A) Pressure at free surface (Figure 17; location point 3). (B) Pressure at bottom surface (Figure 17, location point 6). (b) Pressure time history of S175 hull at stern region, Fn=0.25, λ/L=1. (A) Pressure at free surface (Figure 17; location point 3). (B) Pressure at bottom surface (Figure 17; location point 6).
(a)(b)Figure 21
(a) Contour plot of pressure profile for a particular instant of time of S175 ship,Fn=0.25, λ/L=1. (b) Contour plot of pressure profile for a particular instant of time of ULC, Fn=0.25, λ/L=1.
(a)(b)


## 5. Conclusions

This 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. The 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.



---

*Source: 4564769-2022-03-08.xml*