Experimental and Numerical Study of the Hydroelastic Response of a River-Sea-Going Container Ship

The hydroelastic behaviour of a river-sea-going ship hull is analysed experimentally and numerically. A segmented ship model connected by a steel backbone is tested in regular waves, and its high-frequency vibrations such as springing and whipping responses are identified. The hydroelastic response of the ship is numerically calculated using a hydroelastic time domain method based on strip theory, which is extended to include an improved model of the slamming load. The slamming forces in the bow section are determined using the Modified Longvinovich Model (MLM) instead of the Von Karman model. The vertical motions and wave-induced loads are calculated and compared with the experimental results. The response amplitude operators of the vertical loads and the high-order harmonics are analysed under different speeds, showing good agreement with the experiments. The slamming loads on the bow section of a river-to-sea ship are predicted utilizing the MLM model and compared with the Arbitrary Lagrangian Eulerian algorithm by LS-DYNA and with the measured results.


Introduction
Ship vibration is an important type of response effect for comfort on board and for ship strength. Wave-induced ship hull vibration can be the result of the transient response to slamming loads [1] or the hydroelastic behaviour of the ship hull [2], referred to as whipping and springing respectively, both of which are essential for the consideration of structural safety and integrity.
Hydroelasticity theory was first proposed and applied in the design stage for vessels and large floating structures by Bishop and Price [3] to evaluate the wave-induced responses, including structural deformation, utilizing strip theory, which improved the traditional method regarding the ship as a rigid body. Hydroelasticity has not been a major concern in the calculations of wave-induced vibrations for conventional ships with lengths typically up to 200 m, with the exception of flexible ships such as ones made of composite materials [4]. The recent increase in the sizes of containerships leading to ships with 300-400 m in length has shown that in that size range, hydroelasticity can become important, leading to springing [5,6]. Thus, the hydroelastic response of ship hull vibration should be taken into consideration in the structural safety and integrity assessment at the design stage [7], especially for ship types that are prone to this type of response.
The river-sea-going ship studied here is an unconventional ship type, with high beam/depth ratios, aimed at sailing along the Yangtze River and also in the East China Sea, where it is subject to panel method, while the structure was determined through a 1D FEM method. Hirdaris et al. [26] determined the structural responses of a bulk carrier utilizing the Timoshenko beam idealization and compared the results obtained with a three-dimensional structural idealization. The comparison of the modal vertical bending moment and shear forces were in good agreement for the first few mode shapes. Lakshmynarayanana and Hirdaris [27] coupled Star-CCM+ and Abaqus. The symmetric dynamic responses were assessed and compared with one-way and two-way coupling simulations.
The rigid body response of the river-sea-going ship studied here was determined in the experiments performed at the China Ship Scientific Research Center in Wuxi, as described by Wang et al. [8]. An experimental programme, conducted to study the slamming force in the bow region of this unconventional type of ship [28], illustrated the importance of accounting for the effect of slamming in the hydroelastic hull response. The same ship is now tested at the Wuhan University of Technology (WUT), focusing on the hydroelastic response. The objective of this paper is to investigate the results of this experimental programme to study the hydroelastic effects on the vertical responses of a river-sea-going ship and to compare them with numerical predictions.
In this paper, experiments of segmented model tests are conducted to investigate the characteristic wave-induced vibrations of the hull girder of a river-to-sea ship. A new type of U-shaped section backbone was designed to reproduce the longitudinal distribution of the ship hull stiffness so as to properly scale the vertical responses of the full-scale ship. The hydroelastic response is simulated using the extended version reported here of the hydroelastic nonlinear time domain method combined with the MLM model, which considers the elevation of the water surface. The ship structure is discretized as a non-uniform continuous Timoshenko beam, and its responses are calculated by the FEM method based on the modal superposition method. The total mass and stiffness matrices are derived from the element matrices, taking into account the deflection and rotation of the beam. The hydrodynamic forces are updated at the exact wetted surface. Moreover, the slamming forces and green water are taken into consideration by the MLM model and the Buchner formulation, respectively, which are described in detail in Sections 3.5 and 3.6. Meanwhile, the slamming pressure calculated by the MLM model is compared with LS-DYNA results and with the measured values from the experiments to validate the extended nonlinear hydroelastic time domain method. Moreover, two different numerical methods for the prediction of wave-induced responses are described, as well as the experimental study. A comparison of the numerical results and measurements is illustrated in order to evaluate the springing responses and their effects on the vertical responses.

Experimental Details
The model tests of the river-sea-going ship were conducted at the towing tank of the Wuhan University of Technology to investigate the wave-induced hull girder vibrations, especially the whipping and springing responses, while the linear rigid body response had already been studied by Wang et al. [8]. A scaled model (1:32) of the river-sea-going ship was constructed with FRP (Fiberglass Reinforced Plastics) that was divided into four segments connected with a U-shaped section backbone in order to reproduce the effect of the deformation modes on the loads. The segmented model had the scaled flexible characteristics of the real ship. The space between each segment was 20 mm to avoid contact between the segments during the movement.

Backbone Design and Measurement
The backbone model was designed to simulate the structural dynamic responses of the prototype, while the inertia and hydrodynamic loads on the segmented ship model were transferred to the backbone during the test. The backbone should reproduce the flexibility characteristics of the real river-to-sea ship. The scaled model hull geometry and sectional backbone were designed according to geometric similarity, kinematic similarity, dynamic similarity, modal similarity and structural dynamics principles. The longitudinal distribution of the weight of the model was also designed to be similar to the prototype. The details of the characteristics of the main dimensions of the ship with its pronounced bow flare are described in Table 1, and the body plan of the ship is shown in Figure 1.   The configuration of the segmented model and backbone is shown in Figure 2. In the stern part, the backbone is inclined in order to fit the hull lines. The shapes and detailed properties of the newly designed U-shaped backbone are shown in Figure 3 and Table 2. The section properties of the cross section are calculated by the ANSYS software. The wave-induced loads are measured by rosette strain gauges placed on the backbone at the three cuts, which were 0.25 PP L , 0.5 PP L and 0.75 PP L , respectively, before the aft perpendicular. The accelerometer was mounted on the centre of gravity to measure the vertical accelerations. The pressure sensors are mounted on the bow flare area.  The configuration of the segmented model and backbone is shown in Figure 2. In the stern part, the backbone is inclined in order to fit the hull lines. The shapes and detailed properties of the newly designed U-shaped backbone are shown in Figure 3 and Table 2. The section properties of the cross section are calculated by the ANSYS software. The wave-induced loads are measured by rosette strain gauges placed on the backbone at the three cuts, which were 0.25L PP , 0.5L PP and 0.75L PP , respectively, before the aft perpendicular. The accelerometer was mounted on the centre of gravity to measure the vertical accelerations. The pressure sensors are mounted on the bow flare area.   The configuration of the segmented model and backbone is shown in Figure 2. In the stern part, the backbone is inclined in order to fit the hull lines. The shapes and detailed properties of the newly designed U-shaped backbone are shown in Figure 3 and Table 2. The section properties of the cross section are calculated by the ANSYS software. The wave-induced loads are measured by rosette strain gauges placed on the backbone at the three cuts, which were 0.25 PP L , 0.5 PP L and 0.75 PP L , respectively, before the aft perpendicular. The accelerometer was mounted on the centre of gravity to measure the vertical accelerations. The pressure sensors are mounted on the bow flare area.

Experimental Programme
The tests were conducted under several cases of regular waves, whose parameters are listed in Table 3. The first three cases were designed to identify the high-frequency springing vibrations, while the last one induced slamming in the bow. The wave elevation was measured by a resistive wave probe mounted in front of the mean position of the bow, and the absolute motion relative to the carriage was obtained by the optical tracking device. The model was towed under the towing carriage, which was fixed at the centre of gravity, while the vertical motion was unconstrained as shown in Figure 4.

Theory
The hydroelastic response of the river-sea-going ship is simulated using an extended version of the hydroelastic nonlinear time domain method based on strip theory that was developed herein. The simulations include the responses to slamming loads represented according to direct pressure integration using the MLM model, which considers the elevation of the water surface. The formulations in Rajendran et al. [14] use the Von Karman model, which regards the water pile up as zero.

Equation of Motion
Considering the nonlinearity of the wetted hull geometry, the hydrostatic and hydrodynamic forces are calculated at each time step. The equation of the motion of the ship can be expressed as follows:

Theory
The hydroelastic response of the river-sea-going ship is simulated using an extended version of the hydroelastic nonlinear time domain method based on strip theory that was developed herein. The simulations include the responses to slamming loads represented according to direct pressure integration using the MLM model, which considers the elevation of the water surface. The formulations in Rajendran et al. [14] use the Von Karman model, which regards the water pile up as zero.

Equation of Motion
Considering the nonlinearity of the wetted hull geometry, the hydrostatic and hydrodynamic forces are calculated at each time step. The equation of the motion of the ship can be expressed as follows: [M] ..
where [M] is the global mass as a positive definite symmetric matrix including the added mass, [B] is the global damping matrix including the damping coefficient matrix and [K] is the global stiffness matrix with the symmetry of a semi-definite matrix. The {u} is the nodal displacement vector, which can be expressed as: According to the modal superposition theory, the distortion of the structure can be shown as the summation of the distortion in the principal modes indicated above, where p r is the principal coordinate, {u r } is the principal mode and [D] is the modal matrix. The motion equation can be described in terms of the principal coordinates as: where [m], [b] and [k] are the generalized mass, damping and stiffness matrices of the dry structure. F k (t) is the generalized hydrodynamic force vector, which can be divided into the radiation force F R , Froude-Krylov force F FK , diffraction force F D , hydrostatic force F S , slamming force F SLAM and green water forces F GW . Those hydrodynamic forces are calculated corresponding to the exact wetted surface area at each time step. Those force terms are explained in Sections 3.2-3.6 below.

Radiation Force
With the advantages of fast simulation, the nonlinear radiation forces are implemented by the time convolution method using the memory function K m kr and infinite frequency added mass [A kr ∞ ], which is updated at the instantaneous wetted surface after each time step. The radiation problem is treated nonlinearly and is solved in the frequency domain in the same way as the real wetted hull geometry using boundary elements. The expression of the nonlinear radiation forces is as follows: Substituting Equations (5) and (4) into (3), the equation of motion in the time domain for the flexible hull can be written as: where N is the maximum number of modes considered in the numerical calculations, m kr and c kr are diagonal matrix terms and b kr is a symmetric matrix. C m kr is the radiation restoration coefficients. Finally, the equation of motion was solved using the finite difference method.
The principal mode vector u r was replaced with its z component w r (x). The memory functions K m kr and radiation restoration coefficients C m kr were calculated from the frequency-dependent hydrodynamic coefficients utilizing: The frequency-dependent added mass A kr (ω) and damping coefficients B kr (ω) are derived from the sectional added mass m(x) and damping coefficients N(x), both calculated for the exact wetted surface by the boundary integral method of Sutulo and Guedes Soares [29]: This has advantages over the multi-parameter conformal mapping [30] used earlier and does not have the problem of the irregular frequencies of the Frank close-fit method [31]. Figure 5 presents the sectional heave coefficients at a speed of 10 knots implemented by the boundary element method.
The frequency-dependent added mass This has advantages over the multi-parameter conformal mapping [30] used earlier and does not have the problem of the irregular frequencies of the Frank close-fit method [31]. Figure 5 presents the sectional heave coefficients at a speed of 10 knots implemented by the boundary element method.

Diffraction Force
Considering the speed of advance, the diffraction force can be expressed as explained by Bishop and Price [3]. Furthermore, the sectional added mass and damping should be updated at each time step according to the exact wetted surface.

Diffraction Force
Considering the speed of advance, the diffraction force can be expressed as explained by Bishop and Price [3]. Furthermore, the sectional added mass and damping should be updated at each time step according to the exact wetted surface. where ω and ω e are the wave frequency and encounter frequency, respectively; z is the vertical distance between the centroid of the underwater section and the frame of reference and ζ a is the incident wave amplitude.

Froude-Krylov and Hydrostatic Force
According to the incident wave potential, the Froude-Krylov forces are obtained by integrating the pressure over the wetted surface at each time step: where C x is the wetted cross section contour under the incident wave elevation and dl is the incremental length along the girth of the cross section, n is the unit vector normal to the wetted surface, ρ is the density of the fluid, and g is the gravity acceleration.
The hydrostatic force is assumed as the pressure between the mean water line and the wave crest, which is determined as: where F S static is the static hydrostatic force and H(y, z) is the sectional coordinate relative to the mean water level.

Slamming Force
To determine whether the slamming occurs, the relative vertical motion between the ship and wave surface needs to be monitored to determine the region of the bow that came out of the water and was subjected to the slamming load.
Due to its highly nonlinear nature, the occurrence of slamming is strongly correlated to the vertical impact velocities relative to the wave elevation while they re-entered the water. The relative vertical displacement ξ r and the velocity of a specific section V r , considering the distortion of the structure, can be given as follows: where ξ 3 is the heave motion, ξ 5 is the pitch and ν is the sailing velocity.
The slamming force occurred when the sections re-enter the water after emerging out of the wave profile illustrated above. It is calculated according to direct pressure integration based on the MLM model.
where h(t) is the prescribed penetration depth under the water, c(t) is the half wetted length on the body surface, f (y) is the body shape, f y = f (y)/∂y is the partial derivative of the sectional geometry and d(t) is the splash-up height.
The geometry of the ship sections represented by a polynomial equation and the half wetted length on the section can be represented as: The slamming load from the MLM method can be obtained from Equation (13) considering the elevation of the water surface in the case that the vertical motion is larger than the draft and the relative velocity is beyond the threshold. However, it does not take into account the variation of the added mass in the longitudinal direction.
The slamming load was developed in this hydroelastic time domain method based on the MLM model considering the effect of the water's elevation instead of the Von Karman model, which neglects the local rising up of the surface upon impact, which leads to an underestimated wetted length. The MLM model used in this formulation is compared with the Von Karman formulation adopted by Rajendran and Guedes Soares [32] for two sections as an example of the difference in prediction. The geometry of the section can be represented as polynomial equations, and the fitting data of these two sections are shown in Figure 6. Thus, the parameter can be utilized in Equation (13) where ( ) h t is the prescribed penetration depth under the water, ( ) c t is the half wetted length on the body surface, ( ) f y is the body shape, ( ) y f f y y = ∂ is the partial derivative of the sectional geometry and ( ) d t is the splash-up height.
The geometry of the ship sections represented by a polynomial equation and the half wetted length on the section can be represented as: The slamming load from the MLM method can be obtained from Equation (13) considering the elevation of the water surface in the case that the vertical motion is larger than the draft and the relative velocity is beyond the threshold. However, it does not take into account the variation of the added mass in the longitudinal direction.
The slamming load was developed in this hydroelastic time domain method based on the MLM model considering the effect of the water's elevation instead of the Von Karman model, which neglects the local rising up of the surface upon impact, which leads to an underestimated wetted length. The MLM model used in this formulation is compared with the Von Karman formulation adopted by Rajendran and Guedes Soares [32] for two sections as an example of the difference in prediction. The geometry of the section can be represented as polynomial equations, and the fitting data of these two sections are shown in Figure 6. Thus, the parameter can be utilized in Equation (13) to estimate the slamming load.  The non-dimensional coefficients of the slamming pressure are compared in Figure 7. The Von Karman model underestimated the slamming effects, utilizing the infinite frequency added mass as a function of submergence relative to the undisturbed free surface. Meanwhile, the time history of the pressure was also different from the results of the MLM, which was obtained from neglecting the water rise up, especially in the initial phase. The neglected large rate of change of the wetted surface is one of the most significant deviation in the non-dimensional coefficients of slamming pressure. The non-dimensional coefficients of the slamming pressure are compared in Figure 7. The Von Karman model underestimated the slamming effects, utilizing the infinite frequency added mass as a function of submergence relative to the undisturbed free surface. Meanwhile, the time history of the pressure was also different from the results of the MLM, which was obtained from neglecting the water rise up, especially in the initial phase. The neglected large rate of change of the wetted surface is one of the most significant deviation in the non-dimensional coefficients of slamming pressure. There are other methods and formulations that have been investigated by Wang and Guedes Soares [33][34][35][36][37] and Wang et al. [28], which consider the effect of the local structure's deformation on the slamming force on the bow, which could be considered in future implementations. The wave and whipping bending moments have been correlated by a probabilistic load combination method [38]. The influence of the correlations was further investigated by a comparison of numerical and There are other methods and formulations that have been investigated by Wang and Guedes Soares [33][34][35][36][37] and Wang et al. [28], which consider the effect of the local structure's deformation on the slamming force on the bow, which could be considered in future implementations. The wave and whipping bending moments have been correlated by a probabilistic load combination method [38]. The influence of the correlations was further investigated by a comparison of numerical and experimental results.

Green Force
The effect of the green water forces on the deck can be calculated on the basis of the momentum method. Thus, the vertical force integrated along the length can be given as [39]: in which the mass of the water on the deck m gw can be derived from the vertical distance between the relative motion and the freeboard, while w is the velocity of the deck.

Structural Response
The disadvantage of utilizing a rigid body model for the wave-induced response is that it does not consider the flexible modes and rotary inertia effects of the structure. In order to account for the coupling effects of the deflections of high-order modes, the main structural characteristics of the hull along the ship length needed to be incorporated in the model.
The response of the flexible hull, including the rotary inertia effect and flexible modes, can be determined by a modal analysis. With the assumption that the ship could be idealized as a non-uniform Timoshenko beam, it is discretized into a number of N continuous beam elements, considering the shear deformation and rotary inertia.
The positive and negative bending moments represent the hogging and sagging moments, respectively. The sectional bending moment is derived from the modal bending moment M r (x), adopting:

Numerical Model
The wave-induced ship motions and nonlinear global loads are assessed by comparing the hydroelastic time domain code partially developed here with the experimental results. To also serve as a reference for the rigid body motion calculations, use is made of the commercial software of Sesam-HydroD, which uses the 3D Rankine panel method for the 3D diffraction and radiation problem, but has no hydroelastic capability.
The mesh of the ship and the free surface in HydroD-Wasim are shown in Figure 8. Approximately 86 panels are generated in the longitudinal hull direction, while there are 34 panels along the transverse direction. The free surface meshes are generated with a radius about 10 times the ship length. The body nonlinearity is accounted for by HydroD-Wasim. Beyond that, the Froude-Krylov and restoring hydrostatic forces are solved up to the instantaneous wetted surface considered in the nonlinear effects. However, the radiation and diffraction forces are linearized and calculated for the mean free surface, which means that the degree of the nonlinearity of Wasim is only partial.
The mesh of the ship and the free surface in HydroD-Wasim are shown in Figure 8. Approximately 86 panels are generated in the longitudinal hull direction, while there are 34 panels along the transverse direction. The free surface meshes are generated with a radius about 10 times the ship length. The body nonlinearity is accounted for by HydroD-Wasim. Beyond that, the Froude-Krylov and restoring hydrostatic forces are solved up to the instantaneous wetted surface considered in the nonlinear effects. However, the radiation and diffraction forces are linearized and calculated for the mean free surface, which means that the degree of the nonlinearity of Wasim is only partial.
. Based on Froude scaling, the bending modulus E is scaled with λ , and the area moment of inertia is scaled as 4 λ , where λ is the Froude scale. Similarly, the shear rigidity kGA was scaled as 3 λ , where k is the shear correction factor, which was taken as 5/6, G is the shear modulus, and A is the cross sectional area of the structure. Figure 9 shows the distribution of the mass (ship's scale) used in the finite element model, in which the structure was assumed to be a 2D non-uniform Timoshenko beam. In the hydroelastic time domain numerical simulations, the ship is modelled with 20 sections. Based on Froude scaling, the bending modulus E is scaled with λ, and the area moment of inertia is scaled as λ 4 , where λ is the Froude scale. Similarly, the shear rigidity kGA was scaled as λ 3 , where k is the shear correction factor, which was taken as 5/6, G is the shear modulus, and A is the cross sectional area of the structure. Figure 9 shows the distribution of the mass (ship's scale) used in the finite element model, in which the structure was assumed to be a 2D non-uniform Timoshenko beam. The exact calculation of the structural damping is complicated; therefore, the damping is assumed to be a linear combination of the mass and stiffness matrices. The critical damping ratio is chosen as 2.7% according to the free decay test results. The deflections and slopes up to third flexible modes for the ship calculated by the modal superposition analysis are presented in Figure 10, and the wet natural frequencies calculated for the flexible model are in good agreement with the measured values as shown in Table 4, even though the numerical results were slightly overestimated. There is good agreement between the measured and the numerical values for the lower modes; however, the disparity increases for higher modes. The wet natural frequency of the two-node vertical bending moment occurred at 6.25 and 6.39 rad/s in the experiment and numerical calculation, respectively.

Mode Experiment Numerical
Mass/L(t/m) The exact calculation of the structural damping is complicated; therefore, the damping is assumed to be a linear combination of the mass and stiffness matrices. The critical damping ratio is chosen as 2.7% according to the free decay test results. The deflections and slopes up to third flexible modes for the ship calculated by the modal superposition analysis are presented in Figure 10, and the wet natural frequencies calculated for the flexible model are in good agreement with the measured values as shown in Table 4, even though the numerical results were slightly overestimated. There is good agreement between the measured and the numerical values for the lower modes; however, the disparity increases for higher modes. The wet natural frequency of the two-node vertical bending moment occurred at 6.25 and 6.39 rad/s in the experiment and numerical calculation, respectively.  Figure 10, and the wet natural frequencies calculated for the flexible model are in good agreement with the measured values as shown in Table 4, even though the numerical results were slightly overestimated. There is good agreement between the measured and the numerical values for the lower modes; however, the disparity increases for higher modes. The wet natural frequency of the two-node vertical bending moment occurred at 6.25 and 6.39 rad/s in the experiment and numerical calculation, respectively.

Results
In this section, the experimental results measured in the wave tank and the numerical results are analysed and compared. The numerical results are represented by the acronyms "Htd" and "HydroD" in the following sections to denote the values numerically calculated by the hydroelastic time domain method and the software of HydroD-Wasim, respectively. The acronym "Exp" shows the experimental results.

RAOs (Response Amplitude Operators) of Vertical Responses
The RAOs of the vertical ship responses are analysed for the validation of the extended hydroelastic time domain method. The numerical and the experimental vertical response RAOs are given in Figure 11 against the encounter frequency. In this section, the experimental results measured in the wave tank and the numerical results are analysed and compared. The numerical results are represented by the acronyms "Htd" and "HydroD" in the following sections to denote the values numerically calculated by the hydroelastic time domain method and the software of HydroD-Wasim, respectively. The acronym "Exp" shows the experimental results.

RAOs (Response Amplitude Operators) of Vertical Responses
The RAOs of the vertical ship responses are analysed for the validation of the extended hydroelastic time domain method. The numerical and the experimental vertical response RAOs are given in Figure 11 against the encounter frequency. The numerical results include the hydroelastic time domain results and the Rankine panel method results. The plots on the left and right illustrate the ship responses with speeds of 10 and 14 knots, which correspond to Froude numbers of 0.144 and 0.202, respectively.
The heave RAOs from the two methods of numerical calculations match reasonably well with the experimental results, especially for the middle frequency range. However, the heave responses at the high frequencies are lower than the measurements for the Froude number of 0.202. As the incident waves become steeper, the discrepancy becomes more noticeable due to stronger nonlinearity. Strip theory predicts lower heave amplitudes than the panel code, except at very low frequencies. The results of the pitch RAOs show a similar trend against the encounter frequency, although both numerical results are a bit lower than the test results. However, with the increase in speed, the RAOs of the pitch were reduced.
The non-dimensional VBM reach its peak value for all these three results when the ratio of the wave length to the ship's perpendicular length is around 1.1 at both speeds. Moreover, with the increase in speed, the RAOs are enhanced as well, which is similar to the heave responses. The HydroD, implementing the Rankine panel method, overestimated the responses remarkably, while the hydroelastic numerical method coincided with the experimental results, even if there was still a slight discrepancy at the higher speed. Therefore, from a practical point of view, the hydroelastic time The heave RAOs from the two methods of numerical calculations match reasonably well with the experimental results, especially for the middle frequency range. However, the heave responses at the high frequencies are lower than the measurements for the Froude number of 0.202. As the incident waves become steeper, the discrepancy becomes more noticeable due to stronger nonlinearity. Strip theory predicts lower heave amplitudes than the panel code, except at very low frequencies.
The results of the pitch RAOs show a similar trend against the encounter frequency, although both numerical results are a bit lower than the test results. However, with the increase in speed, the RAOs of the pitch were reduced.
The non-dimensional VBM reach its peak value for all these three results when the ratio of the wave length to the ship's perpendicular length is around 1.1 at both speeds. Moreover, with the increase in speed, the RAOs are enhanced as well, which is similar to the heave responses. The HydroD, implementing the Rankine panel method, overestimated the responses remarkably, while the hydroelastic numerical method coincided with the experimental results, even if there was still a slight discrepancy at the higher speed. Therefore, from a practical point of view, the hydroelastic time domain method provides a compromise between efficiency and accuracy, which is essential, especially during the design stage.

Comparison of Time and Frequency Domain Responses
The time series and spectral density of the measured data and numerical values at each speed are compared in Figure 12 to Figure 13 separately, and they present satisfactory agreement. Both the heave and pitch values are derived at the centre of mass, while the vertical bending moment corresponds to amidships. domain method provides a compromise between efficiency and accuracy, which is essential, especially during the design stage.

Comparison of Time and Frequency Domain Responses
The time series and spectral density of the measured data and numerical values at each speed are compared in Figures 12 to 13 separately, and they present satisfactory agreement. Both the heave and pitch values are derived at the centre of mass, while the vertical bending moment corresponds to amidships.
The heave values obtained from HydroD are greater than the hydroelastic simulations and the measurements in these two regular waves at both speeds. It is obvious that the difference in heave for the wave-to-ship length of 0.9 is higher than the other, from the comparison shown in Figure 11. At same time, pitch shows more consistency among these three results. Meanwhile, with the increase in the encounter frequency resulting from the speed effect, the amplitudes of the motions are also increased, and they show a slight nonlinearity. The discrepancies among these three results are even higher in Case 3 compared to Case 2.   domain method provides a compromise between efficiency and accuracy, which is essential, especially during the design stage.

Comparison of Time and Frequency Domain Responses
The time series and spectral density of the measured data and numerical values at each speed are compared in Figures 12 to 13 separately, and they present satisfactory agreement. Both the heave and pitch values are derived at the centre of mass, while the vertical bending moment corresponds to amidships.
The heave values obtained from HydroD are greater than the hydroelastic simulations and the measurements in these two regular waves at both speeds. It is obvious that the difference in heave for the wave-to-ship length of 0.9 is higher than the other, from the comparison shown in Figure 11. At same time, pitch shows more consistency among these three results. Meanwhile, with the increase in the encounter frequency resulting from the speed effect, the amplitudes of the motions are also increased, and they show a slight nonlinearity. The discrepancies among these three results are even higher in Case 3 compared to Case 2.   The heave values obtained from HydroD are greater than the hydroelastic simulations and the measurements in these two regular waves at both speeds. It is obvious that the difference in heave for the wave-to-ship length of 0.9 is higher than the other, from the comparison shown in Figure 11. At same time, pitch shows more consistency among these three results. Meanwhile, with the increase in the encounter frequency resulting from the speed effect, the amplitudes of the motions are also increased, and they show a slight nonlinearity. The discrepancies among these three results are even higher in Case 3 compared to Case 2.

High-Order Harmonics of Ship Responses
The comparisons between the high-order harmonics of the numerical and experimental VBM up to the seventh order with Froude numbers of 0.144 and 0.202 are shown in Figure 14; Figure 15

High-Order Harmonics of Ship Responses
The comparisons between the high-order harmonics of the numerical and experimental VBM up to the seventh order with Froude numbers of 0.144 and 0.202 are shown in Figure 14; Figure 15, respectively. The black lines show the nonlinear hydroelastic time domain calculations, while the blue lines are extracted from HydroD. Figure 14. Non-dimensional higher harmonics of the vertical bending moments up to 7th order at amidships in regular wave (Fr = 0.144). The plots (a-f) show the second to 7th harmonics, M5 (2) to M5 (7) .
The first order of the harmonics under several encounter frequencies has already been illustrated in Figure 11. The high-order harmonics are normally around 10% of the first harmonics, especially at a low-encounter frequency, avoiding the vicinity where the wet natural frequencies dominate. It should be noted that normally, the harmonics around the natural wet frequency of the model rise dramatically due to the amplifying effect of the springing vibration, which amounts to 31% of the first ones. The nonlinear hydroelastic time domain method is able to simulate the springing with such a strong nonlinear effect. However, even though it captured the high-order harmonics' characteristics, the amplified value still underestimated the test results. Figure 15. Non-dimensional higher harmonics of the vertical bending moments up to 7th order at amidships in regular wave (Fr = 0.202). The plots (a-f) show the second to 7th harmonics, M5 (2) to M5 (7) .
There is a large discrepancy in the sixth and seventh harmonics only when the encounter frequency is around the natural frequency of the model, but away from that frequency, the agreement is fine. The discrepancy being only at the natural frequency reflects some difficulty in the estimation of the corresponding damping, which is always very difficult for those high-order harmonics.
Along with the increase in speed, the second and third harmonics present similar trends with enlarged values. Again, the Rankine panel method could not identify the harmonics higher than the fourth order. The amplified peaks in the vicinity of the wetted natural frequency still show a distinct The first order of the harmonics under several encounter frequencies has already been illustrated in Figure 11. The high-order harmonics are normally around 10% of the first harmonics, especially at a low-encounter frequency, avoiding the vicinity where the wet natural frequencies dominate. It should be noted that normally, the harmonics around the natural wet frequency of the model rise dramatically due to the amplifying effect of the springing vibration, which amounts to 31% of the first ones. The nonlinear hydroelastic time domain method is able to simulate the springing with such a strong nonlinear effect. However, even though it captured the high-order harmonics' characteristics, the amplified value still underestimated the test results.
There is a large discrepancy in the sixth and seventh harmonics only when the encounter frequency is around the natural frequency of the model, but away from that frequency, the agreement is fine. The discrepancy being only at the natural frequency reflects some difficulty in the estimation of the corresponding damping, which is always very difficult for those high-order harmonics.
Along with the increase in speed, the second and third harmonics present similar trends with enlarged values. Again, the Rankine panel method could not identify the harmonics higher than the fourth order. The amplified peaks in the vicinity of the wetted natural frequency still show a distinct nonlinear in the high-order vibrations caused by the springing effects. Moreover, such effects are enhanced with the increase in the speed, while the wave heights do not change remarkably. The differences in the high-order harmonics of these two simulations are caused by several aspects.
Although HydroD-Wasim is a three-dimensional method, its nonlinearity is still simplified. The body nonlinearity is accounted for, the quadratic pressure is considered and the Froude-Krylov force and restoring forces are calculated up to the instantaneous wetted surfaces, but the radiation and diffraction forces are still linear and are solved on the mean free surface and mean wetted surface. Furthermore, the flexibility of the ship structure is not modelled, as it only takes into consideration the mass distribution and not the stiffness characteristics.
On the other hand, the nonlinear hydroelastic time domain code included the effect of slamming and green water forces as nonlinear loads. The nonlinear radiation force could account for the changes in the hull instantaneously wetted surface. The effect of the flexibility of the structure and the hydroelasticity on the ship responses is also modelled.

Slamming Loads and Their Occurrence
Since the slamming load depends very much on the deadrise angle, it is important to assess the slamming load for the river-to-sea ship's cross sections with small deadrise angles. According to Section 3.5, the specific region in the bow where the slamming load occurs could be assessed. For instance, the relative motion and velocity of the location of P1 on Sta.19 are shown in Figure 16, as well as the wave elevation, while the cross section of Sta.19 and the slamming pressure location are presented in Figure 17. nonlinear in the high-order vibrations caused by the springing effects. Moreover, such effects are enhanced with the increase in the speed, while the wave heights do not change remarkably. The differences in the high-order harmonics of these two simulations are caused by several aspects. Although HydroD-Wasim is a three-dimensional method, its nonlinearity is still simplified. The body nonlinearity is accounted for, the quadratic pressure is considered and the Froude-Krylov force and restoring forces are calculated up to the instantaneous wetted surfaces, but the radiation and diffraction forces are still linear and are solved on the mean free surface and mean wetted surface. Furthermore, the flexibility of the ship structure is not modelled, as it only takes into consideration the mass distribution and not the stiffness characteristics.
On the other hand, the nonlinear hydroelastic time domain code included the effect of slamming and green water forces as nonlinear loads. The nonlinear radiation force could account for the changes in the hull instantaneously wetted surface. The effect of the flexibility of the structure and the hydroelasticity on the ship responses is also modelled.

Slamming Loads and Their Occurrence
Since the slamming load depends very much on the deadrise angle, it is important to assess the slamming load for the river-to-sea ship's cross sections with small deadrise angles. According to Section 3.5, the specific region in the bow where the slamming load occurs could be assessed. For instance, the relative motion and velocity of the location of P1 on Sta.19 are shown in Figure 16, as well as the wave elevation, while the cross section of Sta.19 and the slamming pressure location are presented in Figure 17.  nonlinear in the high-order vibrations caused by the springing effects. Moreover, such effects are enhanced with the increase in the speed, while the wave heights do not change remarkably. The differences in the high-order harmonics of these two simulations are caused by several aspects. Although HydroD-Wasim is a three-dimensional method, its nonlinearity is still simplified. The body nonlinearity is accounted for, the quadratic pressure is considered and the Froude-Krylov force and restoring forces are calculated up to the instantaneous wetted surfaces, but the radiation and diffraction forces are still linear and are solved on the mean free surface and mean wetted surface. Furthermore, the flexibility of the ship structure is not modelled, as it only takes into consideration the mass distribution and not the stiffness characteristics.
On the other hand, the nonlinear hydroelastic time domain code included the effect of slamming and green water forces as nonlinear loads. The nonlinear radiation force could account for the changes in the hull instantaneously wetted surface. The effect of the flexibility of the structure and the hydroelasticity on the ship responses is also modelled.

Slamming Loads and Their Occurrence
Since the slamming load depends very much on the deadrise angle, it is important to assess the slamming load for the river-to-sea ship's cross sections with small deadrise angles. According to Section 3.5, the specific region in the bow where the slamming load occurs could be assessed. For instance, the relative motion and velocity of the location of P1 on Sta.19 are shown in Figure 16, as well as the wave elevation, while the cross section of Sta.19 and the slamming pressure location are presented in Figure 17.  The results are obtained by the hydroelastic time domain numerical simulation and by HydoD and compared with the experimental results, even though the HydroD had not taken into consideration the slamming loads. The estimated entry and exit points are marked in the curve of P1's relative motion, with the λ around the 1.1L PP . Meanwhile, the peak velocity is also marked in the curve, which shows that the velocity reached its peak value just after the corresponding location of entry to the water.
It is shown that the measured impact velocities were larger than the prediction results. The calculated motion and velocity relative to the wave surface agree well with the measured results, except for a slight discrepancy for the peak and phase angle. The numerical results for the vertical velocity are lower than the measured ones. The simplifications used in the numerical method are one reason for the discrepancy in the relative motion. Beyond that, the wave elevation at the bow was transferred spatially according to the value at amidships utilizing linear dispersion. The bow slamming phenomena were recorded and are shown in Figure 18. During the slamming process, the bow of the ship model emerged from the water and impacted the water at a relatively high velocity. Under the influence of the slamming load, the nonlinearity of the motion and bending moment was more prominent. The flexible vibration was influenced by the slamming force, which in turn affected the hydrodynamic loads.  The results are obtained by the hydroelastic time domain numerical simulation and by HydoD and compared with the experimental results, even though the HydroD had not taken into consideration the slamming loads. The estimated entry and exit points are marked in the curve of P1's relative motion, with the λ around the 1.1 PP L . Meanwhile, the peak velocity is also marked in the curve, which shows that the velocity reached its peak value just after the corresponding location of entry to the water.
It is shown that the measured impact velocities were larger than the prediction results. The calculated motion and velocity relative to the wave surface agree well with the measured results, except for a slight discrepancy for the peak and phase angle. The numerical results for the vertical velocity are lower than the measured ones. The simplifications used in the numerical method are one reason for the discrepancy in the relative motion. Beyond that, the wave elevation at the bow was transferred spatially according to the value at amidships utilizing linear dispersion. The bow slamming phenomena were recorded and are shown in Figure 18. During the slamming process, the bow of the ship model emerged from the water and impacted the water at a relatively high velocity. Under the influence of the slamming load, the nonlinearity of the motion and bending moment was more prominent. The flexible vibration was influenced by the slamming force, which in turn affected the hydrodynamic loads. The measured slamming pressures on P1 were compared with the numerical and analytical calculation results shown in Figure 19. In the analytical calculation, the slamming pressure was obtained by the MLM model illustrated in Section 3.5. The numerical simulation was conducted with LS-DYNA, in which the impact velocity was the relative vertical velocity in the model tests at the The measured slamming pressures on P1 were compared with the numerical and analytical calculation results shown in Figure 19. In the analytical calculation, the slamming pressure was obtained by the MLM model illustrated in Section 3.5. The numerical simulation was conducted with LS-DYNA, in which the impact velocity was the relative vertical velocity in the model tests at the same instants. The numerical modelling of the slamming problem was conducted by an Arbitrary-Lagrangian Eulerian (ALE) algorithm. The slamming event was simulated as a ship section impacting calm water. For an ALE solver, the fluid is solved using a Eulerian formulation, while the structure is discretized by a Lagrangian approach. A penalty coupling algorithm enables the interaction between the body and the fluids. The remap step in the ALE algorithm applies a Half-Index-Shift advection algorithm to update the fluid velocity and history variables. The interface between the solid structure and the fluids is captured by the Volume of Fluid method. relatively small. The time history of the slamming pressure shows that the "momentum slamming", as the later stage impact, is essential for this ship type. The results calculated by the MLM model for the slamming pressure proved to be reasonable. The differences between the test results and the numerical and analytical calculations are mainly due to the three-dimensional effects in the tests. Beyond that, the variation of the impact velocity and impact of the local angle between the incoming wave surface and ship section also matters in the pressure history for such a blunt body.
It is worth stressing at this stage that the analysis presented here modelled the hydroelastic behaviour of the ship hull but slamming forces were calculated on the assumption of locally rigid structure at the location where slamming occurs. Other formulations [40][41] have analysed the local hydroelastic response, which may be worth considering in further studies. Figure 19. Results for wave elevation, relative vertical motion and corresponding impact velocities of P1's location in the condition of Case 4.

Conclusions
In this paper, the hydroelastic response of a river-sea-going ship is analysed experimentally and numerically. A segmented ship model with a U-shaped cross section backbone is tested to reproduce the wave-induced loads and motions. A hydroelastic time domain method is extended to include the Modified Longcinovich Model (MLM) for calculating the slamming loads in such a blunt ship section. The slamming load on the bow represented by the MLM model is compared with the LS-DYNA simulation and with measured values.
The RAOs of the vertical responses in regular waves shows a good agreement with the extended code, while the Rankine panel method for a rigid hull overestimated the first harmonics slightly. Slight deviation is observed for a range of frequencies where large relative motions occur. The MLM model results and numerical pressures are higher than the measured values. It is observed that the measured pressure includes an initial transient peak value and a relatively long duration with a lower amplitude of "momentum slamming". The numerical results contain only the initial peak duration, while the MLM model results underestimate the "momentum slamming". Due to the high beam/depth ratios of the geometry of the river-sea-going ship, the deadrise angle was relatively small. The time history of the slamming pressure shows that the "momentum slamming", as the later stage impact, is essential for this ship type. The results calculated by the MLM model for the slamming pressure proved to be reasonable. The differences between the test results and the numerical and analytical calculations are mainly due to the three-dimensional effects in the tests. Beyond that, the variation of the impact velocity and impact of the local angle between the incoming wave surface and ship section also matters in the pressure history for such a blunt body.
It is worth stressing at this stage that the analysis presented here modelled the hydroelastic behaviour of the ship hull but slamming forces were calculated on the assumption of locally rigid structure at the location where slamming occurs. Other formulations [40,41] have analysed the local hydroelastic response, which may be worth considering in further studies.

Conclusions
In this paper, the hydroelastic response of a river-sea-going ship is analysed experimentally and numerically. A segmented ship model with a U-shaped cross section backbone is tested to reproduce the wave-induced loads and motions. A hydroelastic time domain method is extended to include the Modified Longcinovich Model (MLM) for calculating the slamming loads in such a blunt ship section. The slamming load on the bow represented by the MLM model is compared with the LS-DYNA simulation and with measured values.
The RAOs of the vertical responses in regular waves shows a good agreement with the extended code, while the Rankine panel method for a rigid hull overestimated the first harmonics slightly. Slight deviation is observed for a range of frequencies where large relative motions occur. Considering the high-frequency vibration effects such as springing, the Rankine panel's results only show the second-order harmonics, while the hydroelastic simulations captures the amplified resonant vibrations consistent with the experiment's high-frequency characteristics. The numerical simulation is able to capture the hydroelastic effect and represent nonlinear springing very well. The influence of springing on the vertical bending moment is significant and not negligible when the encounter frequency is in the vicinity of the first flexible wet natural frequency.
It necessary to pay attention to such high beam/depth ratios of the geometry of the river-sea-going ship when evaluating the slamming pressure. The slamming load results for the bow section of such a ship type calculated by the MLM model proved to be reasonable, and they show that the "momentum slamming", at a later stage after the impact, is essential for this ship type.
Author Contributions: Conceptualization and writing-original manuscript, Y.W., W.W. and C.G.S. All authors have read and agreed to the published version of the manuscript.