Investigation on Ship Hydroelastic Vibrational Responses in Waves

The hydroelastic vibrational responses of a large ship sailing in regular and irregular head waves are investigated numerically and experimentally. A 3D time-domain nonlinear hydroelastic mathematical model is established in which the hydrostatic restoring forces and incident wave forces are calculated on the instantaneous wetted hull surface to consider nonlinear effects of the rigid motion and elastic deformation of the hull in harsh waves. Radiation and diffraction wave forces are computed on the mean wetted surface based on the 3D frequency-domain potential flow theory. The slamming loads are calculated by momentum theory and integrated into the hydrodynamic forces. The 1D Timoshenko beam theory is adopted to model the vibrational structural response and is fully coupled with the presented hydrodynamic theory in time-domain to generate the hydroelastic equation of motion. Moreover, self-propelled segmented model tests were conducted in a laboratory wave tank to experimentally investigate the hydroelastic responses of a target ship in regular and irregular head seas. The numerical and experimental results are systemically compared and analyzed, and the established hydroelastic analysis model turns out to be reliable and effective in the prediction of ship hydroelastic responses in waves.


Introduction
Structural evaluation of marine vessels during their design stage relies on accurate prediction of hydrodynamic loading and subsequent wave-structure interactions. For the cases where the elastic deformation of the structure plays a negligible role in determining the hydrodynamic forces, the structural analysis can be performed separately from the hydrodynamic analysis. However, large vessels with pronounced bow flare are susceptible to frequent and severe wave impact, subjecting them to subsequent significant transient hull vibrations, called whipping, which inevitably interferes with the surrounding fluid field and violates the basic assumption of the conventional decoupled method mentioned above. Hydroelasticity theory sufficiently considers the effects of structural vibrational responses on the hydrodynamic forces, and therefore, provides a framework in which structural and hydrodynamic analyses are fully coupled.
Research on robust and efficient numerical methods for evaluation of ship wave loads and structural responses is still in progress. Bishop and Price [1] first introduced the methodology of hydroelasticity into ship engineering, paving the way for the following fruitful development.
Early research efforts mainly focused on 2D theories in which 2D strip theories, linear or nonlinear, are adopted to analyze the hydrodynamics, while the elastic responses of the ship are represented by the modal superposition method [2][3][4]. With the advent and popularization of high-performance computers, 3D hydroelastic analysis gradually captured the attention of the marine community, especially for the cases where the assumptions of the strip method are not appropriately satisfied [5,6]. The fluid flow is modeled by the 3D potential flow theory, and the ship structure is considered using a 1D or 3D finite element (FE) model [7][8][9]. Recently, Computational Fluid Dynamics (CFD) simulation of fluid flow has been incorporated into the methodology of hydroelasticity [10,11], but the immense amount of time required at present hinders the development of this process, and makes it impractical for use in practical engineering applications.
Together with the development of a theoretical approach is the experimental investigation. Small scale ship model tests in towing tanks have been widely conducted to study hydroelastic responses, such as vertical motions and bending moments, with the experimental conditions ranging from moderate seas to high seas [12][13][14]. Besides those small scale model tests, large scale model tests and even the full-scale ship tests have been conducted to further enhance the understanding of fluid-structure interaction mechanisms [15][16][17][18]. A benchmark study on the performance of seventeen hydroelastic analysis codes applied to a target containership was carried out [19]. The hydroelastic model experiment of the target containership was conducted to further test each seakeeping code. Comprehensive comparisons showed that the participating codes gave scattered results for severe sea conditions, indicating the need for further theoretical research and subsequent verification by experiments. Clauss and Klein [20] experimentally investigated the characteristics of vertical bending moments (M BV ) for different hulls in extreme seas. The analysis showed that large bow flares significantly affect the M BV s.
The understanding and simulation of nonlinearities associated with hydroelastic vibrational responses of large vessels in high seas are the central task for the ship engineering community. The fierce fluid structure interaction is characterized by strongly nonlinear slamming, wave breaking, large amplitude ship motion, etc., posing great technical difficulties on the theoretical prediction of ship responses. More often than not, severe seas have highly-irregular wave conditions which make the characteristics of ship hydroelasticity very different from those in regular waves. Although there has been significant progress in the field of fluid structure interaction, the study of nonlinear hydroelastic responses of large ships in severe irregular sea conditions using 3D methodology is rather scarce in the present literature. Accordingly, the aim of the present paper is to systematically study the hydroelastic vibrational characteristics of a target large ship in both regular and irregular waves. A fully coupled 3D nonlinear hydroelastic method in the time domain is constructed to study whipping responses of a large ship with a pronounced bow flare in regular and irregular wave conditions. The structure of the vessel is idealized as a 1D Timoshenko beam, and the fluid flow is simulated by the 3D potential flow theory. Slamming loads are calculated based on the momentum method. An in-house hydroelastic analysis code is accordingly developed to predict ship responses in regular and irregular waves. To validate the numerical results, experimental tests of a target large ship in regular and irregular head waves were carried out in a towing tank. Quantitative time domain comparisons of numerical and experimental ship responses in regular head seas are discussed, with qualitative time domain comparisons of ship responses, statistical parameters, and probability of exceedance of response extrema being analyzed for ship short-term responses in irregular seas.

Basic Assumption and Description
The general description of ship hydroelastic responses in regular or irregular waves is demonstrated in Figure 1. In the figure, the fluid field Ω F is assumed to be infinite in the horizontal direction, with parameters R, H, S ∞ , and S F respectively being denoted as the infinite radius of Ω F , finite or infinite depth of ΩF, the infinite boundary of ΩF in the horizontal direction, and the free surface of ΩF. SB denotes the wet hull surface. The earth fixed Cartesian frame Ox0y0z0, the translational Cartesian frame oxyz, and the ship hull fixed Cartesian frame Gxbybzb are adopted to study the fluid structure interaction problem. Initially, reference frames Ox0y0z0 and oxyz are spatially identical. Coordinate planes Ox0y0 and oxy are at the calm waterplane with their origins being positioned at the After Perpendicular (A.P.). The translational Cartesian frame oxyz translates at the constant ship velocity U and keeps parallel to Ox0y0z0. The hull fixed reference frame Gxbybzb moves identically with the ship, and G is the Center of Gravity (CG) of the ship. The Galilean transformation relation between these three reference frames at any time t is shown below: where β denotes the azimuth of the ship.

Hydroelastic Equation of Motion
Based on the modal superposition formulation [2], the hydroelastic equation of motion of a large vessel advancing in regular or irregular waves can be expressed as follows:  [FSL] are the Froude-Krylov force matrix, diffraction force matrix, radiation force matrix, hydrostatic restoring force matrix, and slamming force matrix, respectively.
Even though the structural part of Equation (2) is based on the linear theory of solid mechanics, we can still obtain the nonlinear structural responses by practically taking into account the nonlinear effects of some of those five generalized force matrices. Those nonlinear effects should appropriately reveal the real physics of the fluid structure interaction whilst balancing this requirement against the consideration of time efficiency. Therefore, [FI], [FH] and [FSL] are calculated on the instantaneous wet surface of the ship to include nonlinearities associated with hull geometry. In contrast, [FD] and [FR] The earth fixed Cartesian frame Ox 0 y 0 z 0 , the translational Cartesian frame oxyz, and the ship hull fixed Cartesian frame Gx b y b z b are adopted to study the fluid structure interaction problem. Initially, reference frames Ox 0 y 0 z 0 and oxyz are spatially identical. Coordinate planes Ox 0 y 0 and oxy are at the calm waterplane with their origins being positioned at the After Perpendicular (A.P.). The translational Cartesian frame oxyz translates at the constant ship velocity U and keeps parallel to Ox 0 y 0 z 0 . The hull fixed reference frame Gx b y b z b moves identically with the ship, and G is the Center of Gravity (CG) of the ship. The Galilean transformation relation between these three reference frames at any time t is shown below: where β denotes the azimuth of the ship.

Hydroelastic Equation of Motion
Based on the modal superposition formulation [2], the hydroelastic equation of motion of a large vessel advancing in regular or irregular waves can be expressed as follows:  (2) is based on the linear theory of solid mechanics, we can still obtain the nonlinear structural responses by practically taking into account the nonlinear effects of some of those five generalized force matrices. Those nonlinear effects should appropriately reveal the real physics of the fluid structure interaction whilst balancing this requirement against the The following paragraphs will elaborate on each of those five generalized force matrices.

Hydrodynamic Theory in Ship Hydroelasticity
The hydrodynamics in the present paper is based on the potential flow theory [21] which assumes the fluid to be continuous, incompressible and inviscid, with the fluid flow being assumed to be irrotational. In steady state, the fluid velocity potential in reference frame oxyz can be multiplicatively decomposed as [22]: where ϕ is the space-dependent part of ∅, i is the imaginary unit, ω e is the wave encounter frequency which is defined as: where ω is the natural wave circular frequency, k 0 is the wave number. The spatial part of ϕ can be further additively decomposed as: where ζ A denotes the incident wave amplitude, ϕ I denotes the incident wave potential of unit wave amplitude, ϕ D denotes the diffraction wave potential of unit wave amplitude, ϕ Rm denotes the radiation wave potential corresponding to unit motion amplitude of the m-th degree of freedom of the vessel, ν m is the motion amplitude of the m-th degree of freedom, m is bigger than 6 with indices m < 6 denoting rigid degrees of freedom and m > 6 denoting elastic degrees of freedom, respectively. The incident wave potential ϕ I can be obtained by the incident wave field of regular waves. To express those velocity potentials in a more convenient and uniform fashion during the solution procedures, in the following contents we denote ϕ I as ϕ 0 , with ϕ D being denoted as ϕ n+1 . ϕ Rm is replaced by ϕ r where r also ranges from 1 to n which is the same as m. Then we can jointly denote ϕ I , ϕ Rm and ϕ D as ϕ r where r actually ranges from 0 to n + 1. The boundary value problem for the determination of ϕ Rm and ϕ D can be expressed as follows [21]: where n is the normal vector and is defined as positive when it points from the wet boundary surface into the ship interior, u r denotes the r-th mode of vertical displacement of the hull. The present study adopts the zero speed Green's function [23] based boundary element method to numerically solve ϕ Rm and ϕ D . For both regular and irregular wave conditions, the nonlinear hydrostatic restoring forces acting on the elastic hull are calculated as follows: where S(t) denotes the instantaneous wetted body surface in waves, ρ is the fluid mass density, p ka is the k-th principal coordinate, w k denotes the k-th mode of vertical displacement of the hull, F g denotes the sectional hull weight, L is the longitudinal length of the ship.

Incident Wave Force
For ships sailing in regular waves, the nonlinear incident wave force acting on the instantaneous wetted elastic hull is expressed as: The irregular waves can be theoretically generated by the superposition of an infinite number of component regular waves. Based on Equation (8), the nonlinear incident wave forces acting on ships in irregular waves can be expressed as follows: where ζ Ai , ω ei and ϕ 0i are respectively the incident wave amplitude, wave encounter frequency and wave potential of the ith component regular wave.

Diffraction Wave Force
In regular waves, the linear diffraction wave force acting on the wetted elastic hull surface in calm water is expressed as: For ships advancing in irregular waves, the diffraction wave force is obtained by integrating the contribution of each component regular wave that has its corresponding encounter wave frequency. The time-domain convolution integral method is used to obtain the diffraction wave force induced by irregular waves: where h r D (t) is the impulse response function of diffraction wave force for r-th mode of motion, and H r D (iω e ) is the frequency-domain response amplitude operator of diffraction wave force in regular waves. For ships in regular waves, the linear radiation force acting on the mean wetted hull surface in calm water can be written as: where [A(t)] and [B(t)] are called hydrodynamic coefficients matrices of added mass and added damping, respectively. The components of [A(t)] and [B(t)] can be derived as follows [24]: For ships in irregular waves, the wave memory effects which result from the persistence of radiated waves generated by the ship's motion and its scattering of incident irregular waves play an import role in the determination of radiation wave force. The radiation wave force in irregular waves can be expressed as [25]: (15) where [A ∞ ] is the infinite frequency added mass, K(τ) is the retardation function matrix. The component K rk (τ) of K(τ) can be calculated according to the Kramers-Kronig relation [26]: where B rk (ω e ) is the hydrodynamic coefficient of added damping already given in Equation (14). Note that based on the assumption that the wave encounter frequency ω e is relatively high and the vessel velocity U is relatively low [27], the operator U ∂ / ∂x in Equations (6), (8)(9)(10) and (14) can be dropped.

Slamming Wave Force
The slamming loads acting on ships in regular or irregular waves can be obtained by the momentum impact theory which assumes that the slamming force is determined by the momentum changing rate of the fluid around the ship. Therefore, the sectional slamming force can be expressed as follows [28]: where m ∞ (x) is the sectional infinite frequency added mass in heave mode, w rel (x,t) is the hull vertical displacement relative to the incident wave elevation which can be expressed as [2]: where w r (x) is the r-th mode of ship vertical displacement, ζ(x,t) is the incident wave elevation. Then the generalized slamming force can be expressed as [29]: where w(x) is vertical normal vector component of n.

Elasticity Theory in Ship Hydroelasticity
For simplicity, the whole ship structure can be treated as a long non-uniform beam, meaning that the beam theory can be used to describe the ship with the beam sectional parameters, such as the second moment of area, being obtained exactly based on the structural arrangement of the real ship section. This long ship hull beam as a whole can be mathematically described by the Euler-Bernoulli beam theory to predict global structural responses of the ship. However, either by adopting the usual finite element method or simplified linear modal superposition method, the whole ship hull beam needs to be further divided into beam segments, some of which may be short, and therefore, do not satisfactorily fulfill the requirement of adopting Euler-Bernoulli beam theory. Therefore, the present research treats the whole ship hull beam as a segmented Timoshenko beam which takes the effects of shear deformation and rotational inertia on the beam into consideration. The free vertical vibration of the i-th segment of the whole Timoshenko beam can be described by the following differential equation [30,31]: where w i is the vertical displacement, E = 210 GPa and G = 80.77 GPa are respectively the elastic modulus and shear modulus which take the values of Q235 steel (since the typical structural components of the real ship across the chosen ship sections are made from Q235 steel), and are assumed to be invariant across the beam segments; I i is the second moment of area which is considered to be constant along the i-th segment, µ i is the mass per unit length, J i is the rotational inertia per unit length, A i is the cross section area, and k is the Timoshenko shear coefficient and for the rectangular section k = 5/6. We assume that the solution of Equation (20) has the following form: where ω 0i is the natural circular frequency of the vibration, ε is the initial phase angle. Substitute Equation (21) into Equation (20), we have: The solution of Equation (22) is: where each C i is a constant, λ i and the coefficients can be expressed as: In addition to the vertical displacement mode w i (x), the rotation angle mode θ i (x), vertical bending moment mode M i (x) and shear force mode N i (x) can be accordingly obtained, and the solution can be written in matrix form as follows: where: The components of matrix [A] of the right hand side of Equation (25) are constants which can be solved by the boundary condition, and the relation between components of [A] and C i in Equation (23) can be given as follows: Since these beam segments are linked together, we can adopt the Transfer Matrix Method to solve the overall vibrational characteristics of the whole Timoshenko beam of the hull. Note that for the oblique wave conditions, the twisting vibration modes are important and should be considered. Since this study concentrates on ship hydroelastic responses in head wave conditions, the 1D Timoshenko beam theory is adequate because only the vertical vibration mode exists. In addition, by adopting modal superposition formulation, the computational costs between above-mentioned two beam theories, namely the Euler-Bernoulli beam theory and the Timoshenko beam theory, are nearly the same, since the difference only lies in the generalized matrices [m], [c] and [k] in Equation (2), all of which can be conveniently obtained for both beam theories. This fact can also be considered as one advantage of adopting Timoshenko beam theory in the modal superposition method.

Solution of Nonlinear Hydroelastic Equation
For the solution of the hydroelastic equation of motion (i.e., Equation (2)), the fourth-order Runge-Kutta algorithm [32] is adopted in the present study. Then based on the modal superposition principle, the time-domain displacement w(x,t), vertical bending moment M(x,t), and shear force V(x,t) can be written as follows: where w r (x), M r (x) and V r (x) are the r-th mode of vertical displacement, bending moment and shear force, respectively.

Hydroelastic Experiment Setup
A self-propelled segmented model was designed and manufactured according to the target full-scale ship with broad bow flare. Table 1 presents the principal particulars of the ship model. Figure 2 shows the overview of the segmented model. The elastic model is divided into seven segments longitudinally, and sectional loads at the six cut sections are measured during experiments. A total of three steel backbones are used to connect the segments. The backbones are also made of Q235 steel, and therefore, have the same material parameters as those of the real ship hull beam. The heights and widths of these three steel backbones are iteratively adjusted to fulfil the requirements of the principle of similarity on the flexible characteristics for the purpose of reproducing the hydroelastic effect of the real ship. The first order natural frequency of vertical vibration of the ship model after the iteration process mentioned above is adjusted to be similar with that of the real ship. Figure 3 presents the arrangement of the three backbones. The three steel beams are fitted with strain gauges of 200 Hz sampling rate to measure the M BV of the model at the cut sections. A propelling system, positioned at the stern region, is designed to reach the expected speed of the model in still water and waves.   The tests were carried out in the ship model towing tank of Harbin Engineering University. The towing tank is 108 m long, 7 m wide, and 3.5 m deep. A towing carriage is used to serve as the platform where the monitoring apparatus such as the data acquisition system are placed. A computer controlled hydraulic driven wave generator, which is currently utilized as flap type, can generate both regular and irregular waves. The waves will be absorbed by a damping beach opposite to the wave maker to prevent the effects of reflected waves on the generated incident waves and the ship model. A wave probe is positioned near the wave generator to monitor and record the generated waves at 50 Hz sampling rate. Table 2 gives the main parameters of the towing tank. Figure 4 shows the towing tank and the wave generator.  The tests were carried out in the ship model towing tank of Harbin Engineering University. The towing tank is 108 m long, 7 m wide, and 3.5 m deep. A towing carriage is used to serve as the platform where the monitoring apparatus such as the data acquisition system are placed. A computer controlled hydraulic driven wave generator, which is currently utilized as flap type, can generate both regular and irregular waves. The waves will be absorbed by a damping beach opposite to the wave maker to prevent the effects of reflected waves on the generated incident waves and the ship model. A wave probe is positioned near the wave generator to monitor and record the generated waves at 50 Hz sampling rate. Table 2 gives the main parameters of the towing tank. Figure 4 shows the towing tank and the wave generator.    Wave-induced heave and pitch motions of the model are monitored by an in-house designed, four-freedom motion recorder apparatus, which is installed on the towing carriage and connected to the ship model via two heave rods (see Figure 2). The two heave rods are jointly used to measure heave and pitch motions of the ship model while confine sway and yaw motions. The heave of the ship model at its CG can be derived by the two measured heave values at designated positions where the heave rods are installed. The pitch of the ship model is monitored by the inclinometer, which is fixed on the bottom of the heave rod. The sampling rate of motion is set to be 50 Hz during the experiment.

Analysis of Numerical and Experimental Results
The experimental results and the numerical results are compared and analyzed in this section.  Table 3. Note that the two regular wave periods presented in Table 3  Wave-induced heave and pitch motions of the model are monitored by an in-house designed, four-freedom motion recorder apparatus, which is installed on the towing carriage and connected to the ship model via two heave rods (see Figure 2). The two heave rods are jointly used to measure heave and pitch motions of the ship model while confine sway and yaw motions. The heave of the ship model at its CG can be derived by the two measured heave values at designated positions where the heave rods are installed. The pitch of the ship model is monitored by the inclinometer, which is fixed on the bottom of the heave rod. The sampling rate of motion is set to be 50 Hz during the experiment.

Analysis of Numerical and Experimental Results
The experimental results and the numerical results are compared and analyzed in this section.  Table 3. Note that the two regular wave periods presented in Table 3 correspond to the relation of λ⁄L = 1.0, where λ and L respectively denote the wave length and the length waterline of the ship. Figure 5 demonstrates the total panels of the whole ship surface and the partial panels under the calm water surface. Figure 6 presents the extracted cross sections along the hull for slamming loads estimation. Table 3. Experimental conditions extracted for analysis. experiments. The International Towing Tank Conference (ITTC) released guidelines [34] to avoid such interference from the tank side walls in head waves. Based on the guidelines, the two regular wave conditions are free from the tank side wall effects. However, the effect of the side walls is unavoidable for irregular wave conditions because the irregular wave inevitably contains those regular wave components which will induce the side wall effect. To minimize the side wall effect, the first irregular wave condition in the present research has a rather large significant wave height, resulting in severe slamming events and subsequent wave breaking and large eddies that will greatly dissipate the energy of scattered waves. For the second irregular wave condition, only very few regular wave components may cause the side wall effect due to the relatively high vessel speed. In addition, the relatively high ship speed will also give rise to frequent slamming events that will significantly reduce the side wall effect by the energy dissipation mechanism. Therefore, the tank side wall effect in the present study is significantly reduced, and the experimental results are trustable.    The calculation of nonlinear [F I ] and [F H ], as introduced in Section 2.3, is based on the instantaneous wetted hull surface, the acquisition of which can be accomplished by two steps, as illustrated in Figure 7. The first step is to obtain the so-called physical instantaneous wetted surface by intersecting the deformed ship hull panel with the incident wave profile; see the upper part of Figure 7. In this step, special care should be given to those panels whose vertices are neither all below nor all above the incident wave profile. For panels of this kind, the two intersection points between the incident wave profile and the panel sides are joined together by a straight line segment. If the remaining panel is neither a triangle nor a rectangle, then it is further subdivided into a triangle panel and a rectangle panel. After this first step, a satisfactory physical instantaneous wetted surface is acquired because the difference between the straight line segment and the wave profile within the panel is very small. Note that in Equation (6), the boundary of the free surface is actually the calm waterplane; therefore, the physical instantaneous wetted surface obtained by the first step should be mapped to ensure that the clam waterplane is the boundary of the transformed instantaneous wetted surface [33]. Therefore, the second step is the mapping of the physical instantaneous wetted surface, as shown in the lower part of Figure 7, and the mapped instantaneous wetted surface will be adopted.

Necessary Hydrodynamic Results
In this section, hydrodynamic coefficients and retardation functions obtained by the hydroelastic code are presented and analyzed. The added masses of A33 and A55, damping coefficients of B33 and B55, and the retardation functions of K33 and K55 of the target ship with forward speeds 5 kn and 18 kn in head regular waves are presented in Figures 8-10. As can be seen from Figure 8, both A33 and A55 vary significantly with respect to the natural circular frequency ω of the incident wave when ω is less than 0.5 rad/s. For ω greater than 1.5 rad/s, these two added mass coefficients change very slightly. Different from A33, ship velocity plays a more vital role in A55 when ω is less than 0.3 rad/s. As velocity increases from 5 kn to 18kn, positions of occurrence of the peak values in B33 and B55 both move leftward a little bit, with the peak value of B33 changing less significantly than that of B55. Similar to A33 and A55, the variations of B33 and B55 with respect to ω are much more evident when ω falls in the interval of 0.1 rad/s to 1.5 rad/s. The general shapes of retardation functions K33 and K55 are similar, and small oscillations can be found in these two curves. Both K33 and K55 decay quickly with respect to time, as can be observed in Figure 10. Before we proceed to present the experimental and numerical results in the following sections, we briefly discuss the effect of the tank side walls. Compared with the real ship sailing in open seas, the ship model may be affected by the reflected waves from both tank side walls during the experiments. The International Towing Tank Conference (ITTC) released guidelines [34] to avoid such interference from the tank side walls in head waves. Based on the guidelines, the two regular wave conditions are free from the tank side wall effects. However, the effect of the side walls is unavoidable for irregular wave conditions because the irregular wave inevitably contains those regular wave components which will induce the side wall effect. To minimize the side wall effect, the first irregular wave condition in the present research has a rather large significant wave height, resulting in severe slamming events and subsequent wave breaking and large eddies that will greatly dissipate the energy of scattered waves. For the second irregular wave condition, only very few regular wave components may cause the side wall effect due to the relatively high vessel speed. In addition, the relatively high ship speed will also give rise to frequent slamming events that will significantly reduce the side wall effect by the energy dissipation mechanism. Therefore, the tank side wall effect in the present study is significantly reduced, and the experimental results are trustable.

Necessary Hydrodynamic Results
In this section, hydrodynamic coefficients and retardation functions obtained by the hydroelastic code are presented and analyzed. The added masses of A 33 and A 55 , damping coefficients of B 33 and B 55 , and the retardation functions of K 33 and K 55 of the target ship with forward speeds 5 kn and A 55 vary significantly with respect to the natural circular frequency ω of the incident wave when ω is less than 0.5 rad/s. For ω greater than 1.5 rad/s, these two added mass coefficients change very slightly. Different from A 33 , ship velocity plays a more vital role in A 55 when ω is less than 0.3 rad/s. As velocity increases from 5 kn to 18kn, positions of occurrence of the peak values in B 33 and B 55 both move leftward a little bit, with the peak value of B 33 changing less significantly than that of B 55 . Similar to A 33 and A 55 , the variations of B 33 and B 55 with respect to ω are much more evident when ω falls in the interval of 0.1 rad/s to 1.5 rad/s. The general shapes of retardation functions K 33 and K 55 are similar, and small oscillations can be found in these two curves. Both K 33 and K 55 decay quickly with respect to time, as can be observed in Figure 10.

Modal Analysis
Based on the Timoshenko beam theory and the transfer matrix method, the first three vertical vibrational elastic modes are numerically obtained; see Figure 11. Based on the conclusions of [35,36], in the present research we do not give the experimental elastic vertical modes. A hammering test of the experimental model was carried out in calm water to validate the numerical result. The natural frequencies can be obtained by Fourier transform of the decayed time series of the transient structural response following the hammering impact. Table 4 presents the comparison of elastic vibrational wetted natural frequencies between calculation and measurement. The discrepancy between the experimental and numerical results for the wet natural frequency of the first order is 3.439%, which satisfied quite well the requirement of the principle of similarity. Although the errors for the second and third orders of wet natural frequencies are larger than 10%, they have little influence on the overall results, because the effect of the higher order deformation modes on the loads is rather small compared with the first order vibrational mode, which can be verified by Figure 12. In Figure 12, the vertical bending moments obtained based on different number of elastic vibrational modes are presented; it can be concluded that the first order deformation mode is most important, and the contributions from orders of elastic modes greater than three can be safely discarded. A similar conclusion for a linear case was also reported in [37].

Modal Analysis
Based on the Timoshenko beam theory and the transfer matrix method, the first three vertical vibrational elastic modes are numerically obtained; see Figure 11. Based on the conclusions of [35,36], in the present research we do not give the experimental elastic vertical modes. A hammering test of the experimental model was carried out in calm water to validate the numerical result. The natural frequencies can be obtained by Fourier transform of the decayed time series of the transient structural response following the hammering impact. Table 4 presents the comparison of elastic vibrational wetted natural frequencies between calculation and measurement. The discrepancy between the experimental and numerical results for the wet natural frequency of the first order is 3.439%, which satisfied quite well the requirement of the principle of similarity. Although the errors for the second and third orders of wet natural frequencies are larger than 10%, they have little influence on the overall results, because the effect of the higher order deformation modes on the loads is rather small compared with the first order vibrational mode, which can be verified by Figure  12. In Figure 12, the vertical bending moments obtained based on different number of elastic vibrational modes are presented; it can be concluded that the first order deformation mode is most important, and the contributions from orders of elastic modes greater than three can be safely discarded. A similar conclusion for a linear case was also reported in [37].

Analysis of Experimental Regular Waves
The experimental regular waves are generated by the flap type wave maker, and a wave probe was positioned near the wave generator to monitor the waves. Figure 13 illustrates the time history of the experimentally-generated regular wave with the theoretical wave height and wave period being 150 mm and 1.66 s at model scale, respectively. A total of 12 periods are shown in Figure 13, and each wave peak (or wave trough) slightly differs from each other due to the inevitable uncertainty associated with the experimental flap type wave maker. The relatively smooth region of the shown time history is extracted to calculate the mean actual wave height and wave period generated, and the comparison of experimental results and expected values are shown in Table 5. The two regular, experimentally-generated wave conditions are in good agreement with the input regular wave parameters on average, with the biggest relative error of the wave height being 1.81%.

Analysis of Experimental Regular Waves
The experimental regular waves are generated by the flap type wave maker, and a wave probe was positioned near the wave generator to monitor the waves. Figure 13 illustrates the time history of the experimentally-generated regular wave with the theoretical wave height and wave period being 150 mm and 1.66 s at model scale, respectively. A total of 12 periods are shown in Figure 13, and each wave peak (or wave trough) slightly differs from each other due to the inevitable uncertainty associated with the experimental flap type wave maker. The relatively smooth region of the shown time history is extracted to calculate the mean actual wave height and wave period generated, and the comparison of experimental results and expected values are shown in Table 5. The two regular, experimentally-generated wave conditions are in good agreement with the input regular wave parameters on average, with the biggest relative error of the wave height being 1.81%.

Analysis of Experimental Regular Waves
The experimental regular waves are generated by the flap type wave maker, and a wave probe was positioned near the wave generator to monitor the waves. Figure 13 illustrates the time history of the experimentally-generated regular wave with the theoretical wave height and wave period being 150 mm and 1.66 s at model scale, respectively. A total of 12 periods are shown in Figure 13, and each wave peak (or wave trough) slightly differs from each other due to the inevitable uncertainty associated with the experimental flap type wave maker. The relatively smooth region of the shown time history is extracted to calculate the mean actual wave height and wave period generated, and the comparison of experimental results and expected values are shown in Table 5. The two regular, experimentally-generated wave conditions are in good agreement with the input regular wave parameters on average, with the biggest relative error of the wave height being 1.81%.

Motions in Regular Waves
The experimental and numerical hydroelastic responses of heave and pitch at CG of the target ship are extracted and analyzed for the two regular wave conditions. Figures 14 and 15 display the comparison of motion responses in the two testing conditions. The experimental heave and pitch responses show some oscillations in both peaks and troughs compared with the numerical results. This is because the phenomenon of forward speed fluctuation of ships in waves can naturally occur in the present self-propelled ship model experiment, while the present numerical method assumes that the ship velocity remains unchanged during the whole simulation process. In addition, since the experimentally-generated wave peaks (or wave troughs) slightly differ from each other as illustrated in Figure 13, the ship model motion responses under such experimental waves will consequently show some difference in their peak values and trough values, as demonstrated in Figures 14 and 15. The numerical results slightly overestimate the heave and pitch responses on average compared with the experimental results, as can be found in Figures 14 and 15, probably due to the omission of dissipation effects associated with fluid viscosity in the current numerical methodology. Figures 14 and 15 show that, on average, the peak values and trough values of heave response are nearly symmetric. In comparison, both the experimental and numerical trough values of pitch response are larger than peak values. This phenomenon results from the joint effect of nonlinearity of the ship hull geometry and the dynamic effect of waves. Here, negative pitch values represent the trim by stern, meaning that in the present case the specific ship hull geometry tends to lift the ship bow region a little bit in waves. The experimental results of both heave and pitch responses include some minor noise due to the friction and vibration of the motion monitor apparatus. In general, the numerical results agree satisfactorily with the experimental results.

Motions in Regular Waves
The experimental and numerical hydroelastic responses of heave and pitch at CG of the target ship are extracted and analyzed for the two regular wave conditions. Figures 14 and 15 display the comparison of motion responses in the two testing conditions. The experimental heave and pitch responses show some oscillations in both peaks and troughs compared with the numerical results. This is because the phenomenon of forward speed fluctuation of ships in waves can naturally occur in the present self-propelled ship model experiment, while the present numerical method assumes that the ship velocity remains unchanged during the whole simulation process. In addition, since the experimentally-generated wave peaks (or wave troughs) slightly differ from each other as illustrated in Figure 13, the ship model motion responses under such experimental waves will consequently show some difference in their peak values and trough values, as demonstrated in Figures 14 and 15. The numerical results slightly overestimate the heave and pitch responses on average compared with the experimental results, as can be found in Figures 14 and 15, probably due to the omission of dissipation effects associated with fluid viscosity in the current numerical methodology. Figure 14a and Figure 15a show that, on average, the peak values and trough values of heave response are nearly symmetric. In comparison, both the experimental and numerical trough values of pitch response are larger than peak values. This phenomenon results from the joint effect of nonlinearity of the ship hull geometry and the dynamic effect of waves. Here, negative pitch values represent the trim by stern, meaning that in the present case the specific ship hull geometry tends to lift the ship bow region a little bit in waves. The experimental results of both heave and pitch responses include some minor noise due to the friction and vibration of the motion monitor apparatus. In general, the numerical results agree satisfactorily with the experimental results. The results of peaks and troughs of motion responses at different wave heights with forward speed 18kn are demonstrated in Figure 16. As the wave height increases, both the peaks and troughs of the heave and pitch responses increase nonlinearly, with the degree of overestimation of numerical results increasing accordingly. Compared with the peaks of motion responses, the nonlinear effect manifested by the troughs of motion responses is more evident. Figure 16 also shows that, with the increase of wave heights, the asymmetry of peaks and troughs in both heave and pitch responses increases. The asymmetry of peaks and troughs in pitch motion becomes increasingly larger than that in heave motion as the wave height continues to increase. The results of peaks and troughs of motion responses at different wave heights with forward speed 18kn are demonstrated in Figure 16. As the wave height increases, both the peaks and troughs of the heave and pitch responses increase nonlinearly, with the degree of overestimation of numerical results increasing accordingly. Compared with the peaks of motion responses, the nonlinear effect manifested by the troughs of motion responses is more evident. Figure 16 also shows that, with the increase of wave heights, the asymmetry of peaks and troughs in both heave and pitch responses increases. The asymmetry of peaks and troughs in pitch motion becomes increasingly larger than that in heave motion as the wave height continues to increase.

MBV in Regular Waves
Together with the extraction and analysis of motion responses, experimental MBVs amidships of the target ship in the two regular wave conditions are also compared and analyzed with numerically obtained MBVs. The total MBVs are decomposed into low frequency and high frequency components corresponding to forced vibration excited by the regular incident waves and slammed vibration due to wave impact, as presented in Figures 17-19. On average, the numerical total MBVs, low frequency MBVs and high frequency slamming MBVs agree satisfactorily with the experimental counterparts. Similar to the motion responses in Figure 14 and Figure 15, the MBVs predicted by the numerical The results of peaks and troughs of motion responses at different wave heights with forward speed 18kn are demonstrated in Figure 16. As the wave height increases, both the peaks and troughs of the heave and pitch responses increase nonlinearly, with the degree of overestimation of numerical results increasing accordingly. Compared with the peaks of motion responses, the nonlinear effect manifested by the troughs of motion responses is more evident. Figure 16 also shows that, with the increase of wave heights, the asymmetry of peaks and troughs in both heave and pitch responses increases. The asymmetry of peaks and troughs in pitch motion becomes increasingly larger than that in heave motion as the wave height continues to increase.

MBV in Regular Waves
Together with the extraction and analysis of motion responses, experimental MBVs amidships of the target ship in the two regular wave conditions are also compared and analyzed with numerically obtained MBVs. The total MBVs are decomposed into low frequency and high frequency components corresponding to forced vibration excited by the regular incident waves and slammed vibration due to wave impact, as presented in Figures 17-19. On average, the numerical total MBVs, low frequency MBVs and high frequency slamming MBVs agree satisfactorily with the experimental counterparts. Similar to the motion responses in Figure 14 and Figure 15

M BV in Regular Waves
Together with the extraction and analysis of motion responses, experimental M BV s amidships of the target ship in the two regular wave conditions are also compared and analyzed with numerically obtained M BV s. The total M BV s are decomposed into low frequency and high frequency components corresponding to forced vibration excited by the regular incident waves and slammed vibration due to wave impact, as presented in Figures 17-19. On average, the numerical total M BV s, low frequency M BV s and high frequency slamming M BV s agree satisfactorily with the experimental counterparts. Similar to the motion responses in Figures 14 and 15, the M BV s predicted by the numerical method are slightly larger than the results obtained experimentally. The ship model velocity fluctuation and the oscillations of the extrema of the generated waves during the experiment also contribute to the oscillations of hogging and sagging values of both the total M BV s and the component loads. Compared with the experimental results, the present numerical method can also produce asymmetrical peaks and troughs of the total and low frequency M BV s responses resulted from the nonlinearity associated with the ship hull geometry. Figure 18 shows that the experimental and numerical low frequency sagging M BV s agree pretty satisfactorily, with some discrepancies being found between numerical and experimental hogging M BV s which are due to the limited nonlinear factors considered in the numerical method. The peak values of numerical high frequency slamming M BV s are also overestimated compared with the experimental results; see Figure 19. method are slightly larger than the results obtained experimentally. The ship model velocity fluctuation and the oscillations of the extrema of the generated waves during the experiment also contribute to the oscillations of hogging and sagging values of both the total MBVs and the component loads. Compared with the experimental results, the present numerical method can also produce asymmetrical peaks and troughs of the total and low frequency MBVs responses resulted from the nonlinearity associated with the ship hull geometry. Figure 18 shows that the experimental and numerical low frequency sagging MBVs agree pretty satisfactorily, with some discrepancies being found between numerical and experimental hogging MBVs which are due to the limited nonlinear factors considered in the numerical method. The peak values of numerical high frequency slamming MBVs are also overestimated compared with the experimental results; see Figure 19.  fluctuation and the oscillations of the extrema of the generated waves during the experiment also contribute to the oscillations of hogging and sagging values of both the total MBVs and the component loads. Compared with the experimental results, the present numerical method can also produce asymmetrical peaks and troughs of the total and low frequency MBVs responses resulted from the nonlinearity associated with the ship hull geometry. Figure 18 shows that the experimental and numerical low frequency sagging MBVs agree pretty satisfactorily, with some discrepancies being found between numerical and experimental hogging MBVs which are due to the limited nonlinear factors considered in the numerical method. The peak values of numerical high frequency slamming MBVs are also overestimated compared with the experimental results; see Figure 19.       Figure 20 shows the relationship of hogging and sagging values of MBVs with respect to wave heights at the ship speed of 18 kn. As can be seen clearly from Figure 20, the total MBVs, low frequency MBVs and high frequency slamming MBVs all vary nonlinearly with the increase of the wave height. Compared with the hogging and sagging values of the total MBVs and low frequency MBVs, the magnitudes of peaks and troughs of the high frequency slamming MBVs are much closer; see Figure 20c. The evident difference between the absolute values of the hogging and sagging MBVs for the total MBVs and low frequency MBVs is mainly due to the nonlinear effect of the change of wet hull geometry on which the wave pressure acts for a relatively long duration in every incident wave period. In comparison, the high frequency slamming MBVs are caused by transient local wave impact, meaning that during the very short duration of whipping, the change of the wetted hull surface is negligible, and therefore, the above-mentioned effects play little role in the asymmetric characteristics of hogging and sagging MBVs of the high frequency MBVs. The omission of some strong nonlinear effects in extreme waves from the present numerical method, such as wave breaking, large eddies, large ship motions, etc., results in the evident discrepancy between the numerical and experimental results in large to extreme waves. Wave breaking and large eddies generated by water viscosity in experiments will dissipate the energy of the ship model system, leading to mitigated ship model responses. The tank side wall effect may also slightly contribute to the difference between the numerical predictions and the experimental results, since this effect is not considered in the present numerical model. On the other hand, large ship motions in extreme waves indicate that not only the [FI] and [FH], but [FD] and [FR] should be calculated based on the 3D instantaneous wetted ship hull, which is currently still a challenge in potential flow approach. The present numerical methodology can reasonably predict the hydroelastic responses in large to extreme waves, even though some discrepancies can be found. Compared with the low frequency MBVs, the hogging and sagging values of the total MBVs and high frequency slamming MBVs change more rapidly as the wave height increases. This is because the transient high frequency impact MBVs are more sensitive to the local change of the wave impact pressure resulted from the increase of incident wave height. Since the high frequency MBVs constitute an important part of the total MBVs, the total MBVs accordingly change rapidly with the increase of wave heights.

The Irregular Waves
The results in the two typical irregular wave conditions are analyzed in this section. The irregular waves are simulated according to the ISSC spectrum using the expected parameters in Table 6. An example of time series recorded by the wave probe during experiment is shown in Figure 21. Then, the experimental wave spectrum can be obtained by correlation function method using the measured time series. A comparison of the wave spectra between the experimental and target results is shown in Figure 22. The comparison indicates that the experimental waves generated by the wave maker are in fairly good agreement with the expected waves. Moreover, as seen from Table 6, the parameters, i.e., the significant wave heights and characteristic wave periods, of generated waves are also found to be close to the expected values. The results in the two typical irregular wave conditions are analyzed in this section. The irregular waves are simulated according to the ISSC spectrum using the expected parameters in Table 6. An example of time series recorded by the wave probe during experiment is shown in Figure 21. Then, the experimental wave spectrum can be obtained by correlation function method using the measured time series. A comparison of the wave spectra between the experimental and target results is shown in Figure 22. The comparison indicates that the experimental waves generated by the wave maker are in fairly good agreement with the expected waves. Moreover, as seen from Table 6, the parameters, i.e., the significant wave heights and characteristic wave periods, of generated waves are also found to be close to the expected values.

Motions in Irregular Waves
The ship responses during a short-term period of 1 h are simulated by the numerical calculation to obtain enough sample data for the subsequent statistical analysis. Since the length of the towing tank is not long enough for the collection of sufficient experimental sample data in a single run of the ship model, a total of 8 runs of the ship model for each irregular wave condition were conducted. Then, the total time series were obtained by connecting the measured time series of each run in the same sea state, thus providing enough experimental data for the following analysis. The time series of heave and pitch motions at CG of the full-scale ship within 1000 s are qualitatively compared from The results in the two typical irregular wave conditions are analyzed in this section. The irregular waves are simulated according to the ISSC spectrum using the expected parameters in Table 6. An example of time series recorded by the wave probe during experiment is shown in Figure 21. Then, the experimental wave spectrum can be obtained by correlation function method using the measured time series. A comparison of the wave spectra between the experimental and target results is shown in Figure 22. The comparison indicates that the experimental waves generated by the wave maker are in fairly good agreement with the expected waves. Moreover, as seen from Table 6, the parameters, i.e., the significant wave heights and characteristic wave periods, of generated waves are also found to be close to the expected values.

Motions in Irregular Waves
The ship responses during a short-term period of 1 h are simulated by the numerical calculation to obtain enough sample data for the subsequent statistical analysis. Since the length of the towing tank is not long enough for the collection of sufficient experimental sample data in a single run of the ship model, a total of 8 runs of the ship model for each irregular wave condition were conducted. Then, the total time series were obtained by connecting the measured time series of each run in the same sea state, thus providing enough experimental data for the following analysis. The time series of heave and pitch motions at CG of the full-scale ship within 1000 s are qualitatively compared from

Motions in Irregular Waves
The ship responses during a short-term period of 1 h are simulated by the numerical calculation to obtain enough sample data for the subsequent statistical analysis. Since the length of the towing tank is not long enough for the collection of sufficient experimental sample data in a single run of the ship model, a total of 8 runs of the ship model for each irregular wave condition were conducted. Then, the total time series were obtained by connecting the measured time series of each run in the same sea state, thus providing enough experimental data for the following analysis. The time series of heave and pitch motions at CG of the full-scale ship within 1000 s are qualitatively compared from . Note that in the present research, the experimental and numerical irregular waves are independently randomly generated according to the same target ISSC wave spectra, and thus, the experimental and numerical details of the component regular waves, namely, the wave amplitude and initial wave phase of each component regular wave, may be different because they are inherently probabilistic. The elevation of the irregular waves is, in fact, generated by assigning a random initial wave phase to each regular wave component during the numerical simulation. Since the time series of experimental and numerical irregular waves differ in fine details, the consequent experimental and numerical time series of hydroelastic motion responses will accordingly not agree in fine detail as demonstrated from Figures 23-26. In comparison, the experimental and numerical results of heave and pitch responses respectively show similarity in trends and magnitudes, qualitatively verifying the validity of the numerical prediction, but still indicating the need for statistical analysis.
independently randomly generated according to the same target ISSC wave spectra, and thus, the experimental and numerical details of the component regular waves, namely, the wave amplitude and initial wave phase of each component regular wave, may be different because they are inherently probabilistic. The elevation of the irregular waves is, in fact, generated by assigning a random initial wave phase to each regular wave component during the numerical simulation. Since the time series of experimental and numerical irregular waves differ in fine details, the consequent experimental and numerical time series of hydroelastic motion responses will accordingly not agree in fine detail as demonstrated from Figures 23-26. In comparison, the experimental and numerical results of heave and pitch responses respectively show similarity in trends and magnitudes, qualitatively verifying the validity of the numerical prediction, but still indicating the need for statistical analysis. independently randomly generated according to the same target ISSC wave spectra, and thus, the experimental and numerical details of the component regular waves, namely, the wave amplitude and initial wave phase of each component regular wave, may be different because they are inherently probabilistic. The elevation of the irregular waves is, in fact, generated by assigning a random initial wave phase to each regular wave component during the numerical simulation. Since the time series of experimental and numerical irregular waves differ in fine details, the consequent experimental and numerical time series of hydroelastic motion responses will accordingly not agree in fine detail as demonstrated from Figures 23-26. In comparison, the experimental and numerical results of heave and pitch responses respectively show similarity in trends and magnitudes, qualitatively verifying the validity of the numerical prediction, but still indicating the need for statistical analysis. The statistical parameters of the hydroelastic motion responses are adopted to describe the ship responses in short-term irregular waves. The comparison of statistical parameters by measurement and calculation is listed in Table 7. Note that for both irregular wave conditions, we present the mean values of heave and pitch responses without calculating the relative errors. This is because the experimental mean values of both heave and pitch motions contain inevitable experimental uncertainties mainly associated with the initial experimental conditions. For each run, we have to choose the initial position of the ship model as the balance position relative to which the heave and pitch responses are experimentally measured. The realistic water surface in the towing tank after each run always contains small ripples that slightly move the ship model from its theoretical balance position used in numerical simulation. It is this experimental error that makes the mean values of experimental heave and pitch motions greatly deviate from those obtained numerically. However, since this experimental error is negligibly small compared with other statistical parameters, it is still reasonable to compare other numerical and experimental statistical parameters. The results indicate that good agreement is achieved between experimental and numerical significant values for peaks and troughs of motion responses, with the difference being less than 16.84%. The relative errors for peak and trough significant values of both heave and pitch motions appear to be larger for higher speed condition, probably due to the insufficient consideration of nonlinearities associated with ship velocity in numerical simulation. The standard deviation also shows satisfactory results between experiment and calculation, with relative error being less than 13.18%.  The statistical parameters of the hydroelastic motion responses are adopted to describe the ship responses in short-term irregular waves. The comparison of statistical parameters by measurement and calculation is listed in Table 7. Note that for both irregular wave conditions, we present the mean values of heave and pitch responses without calculating the relative errors. This is because the experimental mean values of both heave and pitch motions contain inevitable experimental uncertainties mainly associated with the initial experimental conditions. For each run, we have to choose the initial position of the ship model as the balance position relative to which the heave and pitch responses are experimentally measured. The realistic water surface in the towing tank after each run always contains small ripples that slightly move the ship model from its theoretical balance position used in numerical simulation. It is this experimental error that makes the mean values of experimental heave and pitch motions greatly deviate from those obtained numerically. However, since this experimental error is negligibly small compared with other statistical parameters, it is still reasonable to compare other numerical and experimental statistical parameters. The results indicate that good agreement is achieved between experimental and numerical significant values for peaks and troughs of motion responses, with the difference being less than 16.84%. The relative errors for peak and trough significant values of both heave and pitch motions appear to be larger for higher speed condition, probably due to the insufficient consideration of nonlinearities associated with ship velocity in numerical simulation. The standard deviation also shows satisfactory results between experiment and calculation, with relative error being less than 13.18%.

M BV in Irregular Waves
The time series of M BV s amidships of the full-scale ship within 1100 s are shown in Figures 27 and 28, where the low-frequency components and high-frequency components are also separated from the total M BV s by Fourier filtering similar to the treatment of M BV s in regular wave conditions. Similar to the last section, here the experimental and numerical results are also qualitatively compared, and similar magnitudes and trends between calculation and measurement can be found from comparisons of M BV s. As can be seen, several severe slamming events occurred during ship sailing in irregular waves because of high speed or large wave amplitudes. The slamming M BV s decay rapidly due to the structural damping if two adjacent slamming events do not occur too closely. In irregular wave testing condition 1, the magnitude of slamming M BV s is almost the same as that of low frequency M BV s due to the large wave amplitude of irregular wave condition 1 on average. In comparison, for irregular wave testing condition 2, the magnitude of slamming M BV s is obviously smaller than that of low frequency M BV s due to the relatively smaller wave amplitude of irregular wave condition 2 compared with that of irregular wave condition 1. For these two irregular wave conditions, the experimental occurrence frequencies of wave impact manifested by the slamming M BV s differ slightly, indicating that both large wave amplitude and high ship speed will increase the chance of slamming events. For the shown time series of these two irregular wave conditions, the numerical and experimental occurrence frequencies of wave slamming may differ from each other, since the numerical irregular waves are independently randomly generated. On average, the experimental and numerical time series of M BV s show similar trends and magnitudes. The time series of MBVs amidships of the full-scale ship within 1100 s are shown in Figure 27 and Figure 28, where the low-frequency components and high-frequency components are also separated from the total MBVs by Fourier filtering similar to the treatment of MBVs in regular wave conditions. Similar to the last section, here the experimental and numerical results are also qualitatively compared, and similar magnitudes and trends between calculation and measurement can be found from comparisons of MBVs. As can be seen, several severe slamming events occurred during ship sailing in irregular waves because of high speed or large wave amplitudes. The slamming MBVs decay rapidly due to the structural damping if two adjacent slamming events do not occur too closely. In irregular wave testing condition 1, the magnitude of slamming MBVs is almost the same as that of low frequency MBVs due to the large wave amplitude of irregular wave condition 1 on average. In comparison, for irregular wave testing condition 2, the magnitude of slamming MBVs is obviously smaller than that of low frequency MBVs due to the relatively smaller wave amplitude of irregular wave condition 2 compared with that of irregular wave condition 1. For these two irregular wave conditions, the experimental occurrence frequencies of wave impact manifested by the slamming MBVs differ slightly, indicating that both large wave amplitude and high ship speed will increase the chance of slamming events. For the shown time series of these two irregular wave conditions, the numerical and experimental occurrence frequencies of wave slamming may differ from each other, since the numerical irregular waves are independently randomly generated. On average, the experimental and numerical time series of MBVs show similar trends and magnitudes.   The comparisons of statistical parameters of MBVs for the two irregular wave conditions are listed in Table 8. Similar to the case of statistical analysis for motion responses, here again we give the experimental and numerical mean values of MBVs without computing the relative errors. The reason for that is the same as that for the motion responses. The choice of initial position as the balance position before each run of the ship model unavoidably includes experimental error that makes the experimental mean values greatly deviate from those of numerical results. Compared with the small initial ripples, both these two irregular wave conditions have much larger wave amplitudes, and therefore, the other experimental statistical parameters are still trustable. The negative mean values of numerical simulation for both two irregular wave conditions indicate that on average the sagging MBVs are greater than hogging MBVs because of the ship hull geometry nonlinearity. This is similar to the regular wave conditions discussed in former sections. Experimental and numerical significant values show that, not only the hogging and sagging moments, but also the asymmetric characteristics of hogging and sagging moments in irregular wave condition 1 are bigger than that in irregular wave condition 2. This implies that the wave amplitude plays a more important role in the determination of MBVs than ship velocity. The results in Table 8 indicate that good agreement is achieved between experimental and numerical peak and trough significant values with the difference being less than 15.58%. The standard deviations show acceptable results between experimental and numerical results with the relative error being less than 20.55%.  The comparisons of statistical parameters of M BV s for the two irregular wave conditions are listed in Table 8. Similar to the case of statistical analysis for motion responses, here again we give the experimental and numerical mean values of M BV s without computing the relative errors. The reason for that is the same as that for the motion responses. The choice of initial position as the balance position before each run of the ship model unavoidably includes experimental error that makes the experimental mean values greatly deviate from those of numerical results. Compared with the small initial ripples, both these two irregular wave conditions have much larger wave amplitudes, and therefore, the other experimental statistical parameters are still trustable. The negative mean values of numerical simulation for both two irregular wave conditions indicate that on average the sagging M BV s are greater than hogging M BV s because of the ship hull geometry nonlinearity. This is similar to the regular wave conditions discussed in former sections. Experimental and numerical significant values show that, not only the hogging and sagging moments, but also the asymmetric characteristics of hogging and sagging moments in irregular wave condition 1 are bigger than that in irregular wave condition 2. This implies that the wave amplitude plays a more important role in the determination of M BV s than ship velocity. The results in Table 8 indicate that good agreement is achieved between experimental and numerical peak and trough significant values with the difference being less than 15.58%. The standard deviations show acceptable results between experimental and numerical results with the relative error being less than 20.55%. It is of great significance to investigate the asymmetric characteristics of M BV s in hogging and sagging bending states of the ship especially in severe seas. The comparisons of probability of exceedance of M BV s amidships are shown in Figure 29. As can be observed, the numerical results show the same tendency as the experimental counterparts for all the cases, although with some degree of departure because of the limitation of sample size and other reasons. The sagging M BV s are usually larger than the hogging M BV s due to the high frequency slamming M BV s, especially at large M BV regime. For example, the sagging M BV s are almost twice the values of the hogging M BV s at the probability of exceedance of 10 −3 for irregular wave condition 1. In addition, the nonlinearity of ship hull geometry also contributes to the asymmetric characteristics of the total M BV s, as has been analyzed in the regular wave conditions.

Probability Distribution of MBV
It is of great significance to investigate the asymmetric characteristics of MBVs in hogging and sagging bending states of the ship especially in severe seas. The comparisons of probability of exceedance of MBVs amidships are shown in Figure 29. As can be observed, the numerical results show the same tendency as the experimental counterparts for all the cases, although with some degree of departure because of the limitation of sample size and other reasons. The sagging MBVs are usually larger than the hogging MBVs due to the high frequency slamming MBVs, especially at large MBV regime. For example, the sagging MBVs are almost twice the values of the hogging MBVs at the probability of exceedance of 10 −3 for irregular wave condition 1. In addition, the nonlinearity of ship hull geometry also contributes to the asymmetric characteristics of the total MBVs, as has been analyzed in the regular wave conditions. To further investigate the relationship between the low frequency MBVs and the total MBVs, the corresponding scatter diagrams are drawn in Figure 30, where the hogging and sagging MBVs for both low frequency MBVs (obtained by Fourier filtering mentioned above) and total MBVs in each wave encountering period are counted and presented. The points in first and third quadrants denote hogging and sagging MBVs, respectively. The diagonal that separates both the first and third quadrants into two equal parts denotes that the low frequency MBVs are equal to the total MBVs, meaning that the high frequency slamming MBVs are zero. As illustrated in Figure 30, the high frequency slamming MBVs have more effects on the sagging MBVs than the hogging counterparts. The nonlinear behavior is more pronounced for experimental results than numerical results, especially for the higher speed condition, which can be explained by the fact that only weakly nonlinear effects are considered in the present hydroelasticity theory. To further investigate the relationship between the low frequency M BV s and the total M BV s, the corresponding scatter diagrams are drawn in Figure 30, where the hogging and sagging M BV s for both low frequency M BV s (obtained by Fourier filtering mentioned above) and total M BV s in each wave encountering period are counted and presented. The points in first and third quadrants denote hogging and sagging M BV s, respectively. The diagonal that separates both the first and third quadrants into two equal parts denotes that the low frequency M BV s are equal to the total M BV s, meaning that the high frequency slamming M BV s are zero. As illustrated in Figure 30, the high frequency slamming M BV s have more effects on the sagging M BV s than the hogging counterparts. The nonlinear behavior is more pronounced for experimental results than numerical results, especially for the higher speed condition, which can be explained by the fact that only weakly nonlinear effects are considered in the present hydroelasticity theory.

Conclusions
The present research proposes a 3D time domain nonlinear hydroelastic method to numerically estimate the vibrational responses of a large ship advancing in regular and irregular waves. A self-propelled segmented ship model was designed and tested in a towing tank to validate the numerical method presented in this study. The main findings from the presented experimental and numerical studies are as follows: (1) The presented time domain 3D nonlinear hydroelastic method is capable of predicting the hydroelastic vibrational responses of motions and vertical bending moments of large ships with pronounced bow flares sailing in both regular and irregular waves. The numerical results of motions of heave and pitch and MBVs are satisfactorily validated by the results of self-propelled segmented model tests. (2) By comparing numerical and experimental results of heave and pitch responses in regular waves, evident oscillations in experimental peaks and troughs of motion responses can be detected due to the effect of velocity change and the uncertainty associated with the generation of experimental regular waves. Small asymmetric characteristics in peaks and troughs can be found in motion responses because of the nonlinearity associated with the hull geometry. Moreover, compared with the low frequency and high frequency components, the total MBVs show significant asymmetric characteristics in peaks and troughs resulting from the contribution of ship hull geometry nonlinearity and the high frequency slamming components. (3) Either large wave amplitude or high ship velocity will increase the occurrence frequencies of slamming events in irregular wave conditions. Compared with ship velocity, the wave amplitude plays a more important role in the determination of MBVs. The slamming phenomena have relatively small influence on heave and pitch motions, while they have significant effects

Conclusions
The present research proposes a 3D time domain nonlinear hydroelastic method to numerically estimate the vibrational responses of a large ship advancing in regular and irregular waves. A self-propelled segmented ship model was designed and tested in a towing tank to validate the numerical method presented in this study. The main findings from the presented experimental and numerical studies are as follows: (1) The presented time domain 3D nonlinear hydroelastic method is capable of predicting the hydroelastic vibrational responses of motions and vertical bending moments of large ships with pronounced bow flares sailing in both regular and irregular waves. The numerical results of motions of heave and pitch and M BV s are satisfactorily validated by the results of self-propelled segmented model tests. (2) By comparing numerical and experimental results of heave and pitch responses in regular waves, evident oscillations in experimental peaks and troughs of motion responses can be detected due to the effect of velocity change and the uncertainty associated with the generation of experimental regular waves. Small asymmetric characteristics in peaks and troughs can be found in motion responses because of the nonlinearity associated with the hull geometry. Moreover, compared with the low frequency and high frequency components, the total M BV s show significant asymmetric characteristics in peaks and troughs resulting from the contribution of ship hull geometry nonlinearity and the high frequency slamming components. (3) Either large wave amplitude or high ship velocity will increase the occurrence frequencies of slamming events in irregular wave conditions. Compared with ship velocity, the wave amplitude plays a more important role in the determination of M BV s. The slamming phenomena have relatively small influence on heave and pitch motions, while they have significant effects on M BV s. The high frequency M BV s caused by slamming largely increase the total M BV s on the basis of low frequency M BV s, contributing greatly to the asymmetric characteristics of the total M BV s in irregular wave conditions. Therefore, the slamming loads should be given more attention at ship design stage for structural design and evaluation. (4) For irregular wave conditions, the current research focuses on the statistical parameters of ship hydroelastic responses, since the experimental and numerical irregular waves are just compared statistically. The future study should compare the numerical and experimental irregular waves in fine details of their time series by monitoring the near-ship experimental irregular waves, and therefore, the experimental and numerical ship vibrational responses can be compared in fine details. In this way, some uncertainties associated with the randomly-generated irregular waves can be reduced, which currently results in the evident difference between experimental and numerical statistical parameters of ship responses. In addition, the evident difference between the numerical and experimental mean values of ship responses in irregular waves indicates that, before each run of the ship model, the water surface should be sufficiently calm to reduce the error associated with the choice of the balance position of the ship model.