Numerical investigation into combined global and local hydroelastic response in a large container ship based on two-way coupled CFD and FEA

In this study, for efficient safety assessment of a ship’s hull girder subjected to a combination of the global vertical bending moment (VBM) and the local double-bottom bending moment (DBM), a numerical method based on a coupled computational fluid dynamics (CFD) and finite-element analysis (FEA) is first developed. The mutual two-way coupled CFD and FEA in which the elastic deformation of the ship is taken over to the CFD solutions is employed. The present two-way coupled method is developed in two steps, in a weakly coupled manner and a strongly coupled manner. The developed two-way coupled methods are validated via a comparative study against the available experimental results and the straightforward one-way coupled CFD and FEA in terms of the rigid body motion, VBM, and local water pressure. Using the present strongly coupled method, the added mass with regards to the elastic deformation of the ship is reflected in the results of a decrease in the frequency of 2-node vibration mode of the ship or the vibratory (hydroelastic) component in the water pressure. Then, the DBM is evaluated using the estimated water pressure on the ship bottom. An investigation is finally made into the effect of the hydroelastic component in the DBM on the total DBM and the hull girder ultimate strength.


Introduction
Assessing the safety of a ship's hull girder is one of the concerns among ship designers. The global VBM is among the main structural loads exerted on the ship's hull girder. Therefore, an accurate assessment of the VBM is of particular importance. For recent container ships, the hydroelastic response due to the slamming impact force, namely whipping response, is superposed to the normal wave-induced force, and then, it may enlarge the VBM in ships. Hence, a quantitative estimation method for predicting the slamming impact force and resulting whipping component in the VBM is highly needed.
On the other hand, concerns about the local bending moment, in particular the DBM, have been raised [1]. The DBM may further endanger the structure in addition to the whipping VBM. Several preceding studies demonstrated that the DBM reduces the structural capacity of the hull girder under the hogging condition to some extent [2,3]. These works indicated that the water pressure acting on the bottom induces the DBM and thereby increases the compressive force on the outer plate of the double-bottom structure. It is apparent that a reasonable consideration of a combination of the normal wave-induced VBM, whipping VBM, and the DBM should be made in assessing the structural safety of recent large container ships. Then, a consistent numerical method for predicting these structural loads needs to be developed.
For practical design stages, the potential theory-based codes such as the linear/nonlinear strip theory or the 3D panel method have been widely utilized for both the hydrodynamic force and hydroelastic response evaluations [4,5]. However, it is often pointed out that such potentialbased methods are hard to calculate the three-dimensional distribution of the slamming impact force, albeit the low computational burdens. According to the notable advance of the computational performance in recent years, the CFD has been attracting attentions as a method for evaluating both slamming impact force and whipping VBM. Precursory efforts for utilizing the CFD were made by Seng [6], who developed a numerical method for the hydrodynamic and hydroelastic evaluation of a ship using OpenFOAM [7]. In this work, the Timoshenko beam model was implemented in the code, with which the whipping VBM evaluation is made possible. Several research efforts to utilize the FEA for whipping evaluation have been made. Ley et al. [8] and Oberhagemann [9] developed one-way and two-way coupled methods between the CFD and FEA where a 1D beam was used as a structural FE model and then adopted them to whipping evaluations. They demonstrated that the coupled CFD-FEA techniques could well predict the combined wave-induced and whipping VBM in both regular and irregular head sea conditions. One of the advantages of adopting the CFD would be that it enables us to estimate threedimensional distribution of the hydrodynamic loads. In this view, some attempts to combine the CFD and FEA, where a 3D ship is used as a structural FE model, were also made by El Moctar et al. [10] and Ley and Moctar [11]. They established the one-way coupled CFD and 3D FEA leveraging acoustic elements in order for the whipping VBM to take into account the added mass with regards to the elastic deformation. Although the realistic natural frequency of the wetted ship hull was observed in the whipping VBM, there may be room for further improvement, since the added mass was not accurately evaluated using the acoustic elements. An efficient added mass evaluation using the acoustic elements may remain an open issue, but it will be attained by establishing the two-way coupled CFD and 3D FEA.
Our preceding efforts were also devoted to making use of the CFD for the whipping evaluations. We combined the CFD and FEA using a 3D ship FE model, and developed a numerical method to estimate the whipping VBM [12,13] in a container ship. Then, we applied the developed method to evaluate the combined VBM and DBM [14]. However, since these previous studies were based on the one-way coupled method where the hydrodynamic loads are calculated premising the rigid body ship, the added mass with regards to the elastic deformation could not be represented. To address this issue, further extensions of the method towards the two-way coupled manner are expected.
Provided that an efficient two-way coupled CFD-FEA method is established, the next concern is an evaluation of the DBM using the method. In our previous studies [14,15], a series of evaluations of the combined VBM and DBM were made by adopting the one-way coupled CFD-FEA. However, further attention should be paid to the effect of the flexibility of the ship, namely hydroelastic effect, on the DBM. Tatsumi and Fujikubo [2] and Kawasaki et al. [16], amongst others, made the hull girder strength assessment considering the DBM, however, without considering this effect. As the DBM can be regarded as a local response to water pressure, a consideration of the hydroelastic effect into the DBM should be achieved by utilizing the water pressure taking into account the elastic deformation of the ship.
In the present paper, the hydroelastic effect of the combined VBM and DBM in a ship is investigated numerically. To this end, the development and validation of the two-way coupled method between the CFD and 3D FEA are conducted first for evaluating both the whipping VBM and the DBM. A recent container ship running over the regular wave in head sea condition is simulated. The two-way coupled CFD-FEA is developed in two stages; in a weakly coupled manner and a strongly coupled manner. These two-way coupled methods are validated through a series of comparative studies in terms of the RBM, VBM, and water pressure by comparing with the available experiment and the oneway coupled method. Then, the DBM is evaluated using the estimated water pressure on the ship bottom. Finally, an investigation is made into the effect of the hydroelastic component in the DBM on the total DBM and hull girder ultimate strength (US).

Subject ship
A POST PANAMAX size containership (6600 TEU) is chosen as a subject ship. A series of the towing tank test using the scaled model was conducted in National Maritime Research Institute (NMRI) in 2010 [17]. A schematic of the model is given in Fig. 1. The scaled model was composed of six segments rigidly fixed to an aluminum backbone which was used for the reproduction of vertical bending flexibility of the ship. The VBMs were measured along the backbone. The principal particulars of the model are described in Table 1, while the properties of the backbone are described in Table 2.
In this study, the head sea condition under the regular wave with the design speed (Fn = 0.179) is targeted where the wave height (H w ) and non-dimensional wave length (λ/L pp ) are 10 m (in full scale) and 1.0, respectively. The natural frequency of 2-node vertical bending vibration was found at 9.2 Hz in dry condition and 6.8 Hz in wet condition with the logarithmic damping ratio of 0.076. These values are obtained from a series of hammering tests for the ship model, cf. Ref. [14]. For full scale, respective natural frequencies of 2-node vertical bending vibration are 0.946 Hz in dry condition and 0.699 Hz in wet condition, based on the similarity rule.

CFD model
The CFD based on the finite volume method (FVM) is used for the hydrodynamic force evaluations exerted on the ship hull. Fundamental methodologies on the CFD calculation are referenced to our previous studies [14]. The commercial solver STAR-CCM+ 10.06.009-R8 [18] is adopted for the CFD computations. The CFD computations are conducted in the model scale, and then, the hydrodynamic pressures over the discretized hull surface meshes are calculated. Hereafter, all the results and relevant parameters concerning the CFD will be presented at the full scale instead of the model scale. Figure 2 provides the solution domain and mesh discretization of the present CFD model. Half-width of the solution domain is set 0.8L pp which may be enough to avoid reflected surface waves from the Y-direction disturbing flows near the ship model [6]. The non-uniform meshes, where relatively finer meshes are used near free surface and hull surface regions, were to simulate impact pressure on the hull with sufficient accuracy. The basic strategy of the present CFD meshing comes from our previous research, see Ref. [14]. The outlet boundary condition, i.e., the pressure is set to the hydrostatic pressure of the wave to prevent backflow, is used for the downstream plane boundary of the aft end (X = 0). Regarding the other five boundary planes, the velocity-inlet (Inlet) boundary condition, where the velocity of the wave is prescribed to avoid gradient generated from the wall and flow, is adopted [19]. The Euler Overlay Method (EOM) [20,21] is applied to the CFD solution to reduce the effects from the reflection of surface waves from the boundaries or the free surface disturbance due to the ship motion. The region of the damping zone in applying EOM is defined as is the case in Ref. [14]. Semi-implicit pressure linked equations (SIMPLE) algorithm [22] is employed to achieve an implicit coupling between pressure and velocity in the CFD stage. In the SIMPLE algorithm, the pressure field is calculated as follows: where p n+1 is the updated pressure, p n is the current pressure, p′ is the cell values of the pressure correction, N is the maximum number of iterations, and α is the under-relaxation factor. The value of α ranges from 0.2 to 0.5 conventionally. In this study, α = 0.3 is adopted for the CFD computations unless otherwise specified. The maximum number of iterations, N, is set 5.
Further description of the CFD model is indicated in Table 3. The mesh resolution of the CFD over the free surface region, Δx (longitudinal direction), Δy (transverse direction), and Δz (vertical direction), has been described in Table 1. Configuration of the time increment adopted in this study is the same as that provided in our previous study [14]. The mesh size on the hull surface is 1.77 m in full scale, which is almost the same size as the present FE model as described later.

FE model
Here, the FE modeling adopted in this paper is explained. All the relevant parameters concerning the FEA will be presented at the full scale, while the FEA computations are conducted at the model scale. The mesh discretization of the ship model for the FEA is shown in Fig. 3. The fundamental modeling strategies may be referenced to those given in Ref.
[14]. The dynamic explicit solver implemented in LS-DYNA R9.1.0 is used for the FEA computations. The hull surface is modeled with 1.9 m × 1.9 m shell elements. The waterproof sheet is not modeled in the present FE model, as its rigidity is negligibly small (very thin plastic), and then will not influence the VBM. A U-shaped backbone representing the material characteristics and moment of inertia of area of the experimental model has been inserted inside the ship hull. 13639 elements and 14063 nodes are used in total.
To suppress surge motion of the FE model, displacements of nodes in longitudinal direction (in Fig. 3, X-direction) in the vicinity of gravity center are constrained. Rotations around X-and Z-directions are constrained over the hull, so that unintentional local torsional mode is suppressed. The natural frequency of 2-node vertical bending vibration is found at 0.894 Hz in dry condition. In performing FEA, the gravity acceleration is applied to the FE model in the whole solution period, since the hydrodynamic forces from the CFD are given considering the gravity acceleration. Rayleigh type damping is introduced in the FE model applying the logarithmic damping ratio of 0.076. Each calculation time step in the FEA is 1.9 × 10 −6 s, which is far smaller than the present CFD calculation.

One-way coupled method
The one-way coupled CFD-FEA is used for the simulation and validation of the wave-induced load including the whipping vibration of the ship, which is previously developed in Ref. [14]. The procedure on the one-way coupled method is schematically explained in Fig. 4. In coupling the CFD and FEA, pressure values on the hull surface derived from the CFD are directly applied to the FE model. As the RBM of the ship is solved in the CFD phase, pressure values from the CFD include the hydrostatic components and the hydrodynamic components with inertia forces from the RBM. Here,  To exert the pressure calculated over the CFD grids on the FE meshes, mapping processes are needed as the mesh discretization on the hull surface differs among CFD and FEA. In this study, inverse distance weighting (IDW) method [23] is used as a mapping algorithm. The IDW algorithm is described as follows: where P(x, t) is a pressure value at target FE node of interest at time t, P(x i , t) are pressure values at CFD grids at time t, d(x, x i ) are each distance between CFD grids and target FE nodes, and N grid is the number of CFD grids. In this study, four grids closest to the target FE node are employed to execute mapping analyses, i.e., N grid = 4. As the present IDW algorithm merely selects nearest four grids to calculate mapped pressure values, the spatial smoothness of the mapped pressure is not necessarily assured. However, the local non-smoothness will not influence the global hull girder bending behavior if the integrated force is transferred accurately. Figure 5 provides a comparison of the integrated vertical force time history acted on the bow part segment, segment 1 (see in Fig. 3), calculated from the CFD and mapped result. As found from Fig. 5, the mapped vertical force using the IDW is almost comparable to the CFD result. Moreover, slight differences on the vertical force are further corrected at each coupling time step in the same manner as Ref. [14], to keep the force equilibrium between CFD and FEA. The force correction is applied to six hull segments Fig. 3. The coupling calculation is initiated after computing 15 wave cycles in the CFD stage, to exclude the transient part.
The comparisons of the RBM, i.e., heave and pitch motions, between the CFD and FEA phases using the present one-way coupled method are shown in Fig. 6. As there are unbalance forces between imported pressure and acceleration due to different mesh discretization in the CFD and FEA, reaching the identity of RBM between CFD and FEA is still difficult, and therefore, the RBM found from the FEA seems to include a drift component. Meanwhile, when the VBM time series calculated from the present one-way coupled method is focused, see Fig. 7, the VBM time series is found to be stable even under the RBM including the drift. From this result, one can judge the VBM results from the one-way simulation are reliable, i.e., from 22.461 s.

Two-way coupled method
In general, the two-way coupled CFD and FEA falls into the 'partitioned approach', as the fluid solver and the structure solver are alternatingly integrated. In the partitioned twoway coupled method, information about the force and displacement is transferred via the hull surface. The most basic method in the partitioned approach is the so-called weakly   coupled method. An elementary but popular procedure for the weakly coupled method is the conventional serial staggered (CSS) procedure [24], where the coupling solution between the fluid and the structure solver is advanced once per coupling time step. Indeed, the weakly coupled method guarantees less computational efforts, but it is also well known that the weakly coupled method occasionally raises instabilities due to the so-called artificial added mass effect [25]. Such problem occurs in particular when the fluid-structure interaction (FSI) problems concerning the incompressible fluid with flexible structures suffered from large deformations are solved. The artificial added mass effect is considered to be originated from a fact that no rigorous equilibrium between the fluid and the structure domain is established in each coupling time step, which may be inevitable when the weakly coupled method is adopted. To overcome this drawback, the strongly coupled method has been devised and validated by several researchers [26,27].
In the strongly coupled method, sub-iterations between the fluid and structure solvers are implemented in each coupling time step to find convergences between the fluid and structure implicitly. Using the strongly coupled method, solution stability and accuracy will be improved over the weakly coupled method at the expense of computational efforts. In this study, the weakly and strongly coupled methods between the CFD and FEA are developed and compared alongside each other. Stand-alone interface programs to achieve data exchange between the CFD and FEA solvers, without accessing the core part of these softwares (STAR-CCM+ and LS-DYNA), are developed. For overall control of coupling simulations, i.e., call mapping programs and control programs for STAR-CCM+ and LS-DYNA, a control program written in python is used. A java macro is used for controlling the STAR-CCM+ computations, i.e., initiate, interrupt, and resume the solver. LS-DYNA execution commands are included in the control program written in python. Mapping programs for pressure or displacement, written in Fortran, are invoked after computing STAR-CCM+ or LS-DYNA at each coupling time.

Weakly coupled method
The weakly coupled CFD-FEA is assembled by combining the CFD and FEA solutions in a staggered manner. The calculation procedure of the weakly coupled method according to the CSS of which cycle can be described as follows is employed [24]: 1. Advance the CFD solution to the next time step (t = t n ) using the updated fluid meshes. 2. Convert the forces from the CFD solution at t = t n to the external pressures on the FE model via a mapping process. 3. Advance the FEA solution to t = t n by applying the updated external pressures to the FE model, taking over the previous FEA solution at t = t n − Δt. 4. Transfer the calculated hull deformation from FEA to the CFD mesh and then update the fluid mesh via a mesh morphing process. 5. Repeat from steps 1 to 4. A flowchart of the weakly coupled scheme is shown in Fig. 8. The mesh morphing solver implemented in STAR-CCM+ is adopted to reflect the elastic deformation of the ship hull derived from the previous FEA solution. The CFD in each coupling time step is computed considering both the local and global velocity, which are taken over from the previous step. Through the data exchanging processes as implied by red arrows in Fig. 8, the elastic deformation estimated from the FEA in previous time step is reflected to the solution of the CFD at current time step.
The FE model used in this study is not constrained anywhere to suppress the RBM, and hence, deformation results from the FEA include both the RBM and the elastic deformation. Ideally, the RBM is to be calculated only in FEA then feeding both the elastic deformation and RBM back to CFD. However, a divergence of the coupling simulation occurred at a very early stage. This is because of that the mesh movement between each coupling time step was too large for the present morphing solver to handle it. This problem in the morphing solver has been also pointed recently, by Lakshmynarayanana and Temarel [19]. On the other hand, it can be assumed that the effect of discrepancy in the RBM on elastic deformation would be small, as found from the one-way coupled results described earlier, if the external pressure is properly calculated using deformed ship in the CFD phases. To overcome this divergence problem, the deformation obtained from the FEA is once decomposed into heave, pitch motion, and the elastic component at each coupling time step, and then, only the elastic deformation is fed into morphing analyses. By taking this means, one can circumvent the significantly large mesh movement between each coupling time step, thereby the stability of morphing analyses can be enriched. A noteworthy work has been recently presented in Ref. [19] to overcome this problem, by employing 'overset' grid technique to the CFD model. Figure 9 illustrates the decomposition method schematically. Increments of the heave motion at each coupling time step are concisely calculated by translational displacement of the center of gravity (CoG) of the model. It means: where Δh denotes the heave motion increment at t n + Δt, m i means the mass of each FE node, and z i means the vertical position of each FE node. Meanwhile, the increment of the pitch motion Δθ is calculated by applying the linear approximation method, i.e., the least-squares method, to the vertical displacement of FE nodes along the backbone. Consequently, the elastic deformation component at t n + Δt is derived from the following: Mapping analyses are carried out after each CFD or FEA solution to convey the force or the elastic deformation to the FEA or CFD input at the next coupling time step. The IDW algorithm, Eq. 2, is exploited in these mapping analyses. For instance, in the case of the elastic deformation feedbacks (from FEA to CFD), P(x, t) and P(x i , t) in Eq. 2 shall be reread as a displacement value at target CFD node and displacement values at FEA nodes in each coupling time step, respectively.
A sigmoid function is adopted for interpolating the pressure values between consecutive coupling time steps to avoid applying abrupt pressure change to the FE model. The kth interpolated pressure field between time t n and t n + Δt can be described as follows: where N in denotes the number of interpolating points and P denotes the pressure field. N in = 10 is adopted in this study. To maintain the stability of the coupling solution, the elastic deformation of the FE model is suppressed at the early stage of the coupling simulation, i.e., up to 2.0 s, by increasing the damping ratio in the model.
Coupling time step size Δt is held to that of the CFD computation, i.e., Δt is equal to 0.019 s in full scale, see Table 2. When the weakly coupled method is computed using the under-relaxation factor α = 0.3 in the CFD phases, divergence of the solution occurred, presumably due to the artificial added mass effect. This divergence problem was resolved by reducing the magnitude of α to 0.05. Thus, the results on the weakly coupled method presented later are those when α = 0.05 is used.

Strongly coupled method
The strongly coupled CFD-FEA is also assembled in this study. The fundamental approach to assemble the strongly coupled method can be referenced to the algorithm presented by Storti et al. [28]. The generic cycle of the strongly coupled scheme can be described as follows: 1. Advance the CFD solution to the next time step (t = t n ) using the prescribed fluid mesh and then obtain the load field at kth sub-iteration (F tn k ). 2. Convert the load field F tn k to the external pressure field (P ex,tn k ) on the FE model via a mapping process. 3. Advance the FEA solution to the next time step (t = t n ) using the external pressure field P ex,tn k and then obtain the displacement field at kth sub-iteration (u tn k ). 4. Transfer the displacement field u tn k to the CFD mesh then update the prescribed fluid mesh (X tn k ) via a mesh morphing process. 5. Check the convergence of F tn k and/or u tn k . If the solution is converged, determine the current fields as F tn = F tn k , X tn = X tn k . If not, update k = k+1 and repeat from step 1 to 5. 6. Advance the CFD solution to the next time step using the determined fluid mesh X tn .
Note that if only one sub-iteration is implemented, the above-mentioned procedure becomes the weakly coupled method. As regards the convergence of the force and displacement fields at step 5, there are several ways to find Rigid body motion those convergences. As the present approach adopts distinct fluid and structure solvers, in this study CFD and FEA, finding the absolute converged fields of the pressure and deformation is extremely difficult. Among the usual strongly coupled approaches, it is general that convergence criteria for the force or displacement field are defined then the sub-iteration stages are repeated until the convergence condition is satisfied [29]. However, it can also be estimated that huge computational efforts might be required in case that the computation does not converge quickly. Thus, in this study, the so-called predictor-corrector method [30] is adopted to reach a convergence at each coupling time step while decreasing the number of sub-iteration processing. To facilitate the convergence, the Aitken's accelerator [31] is introduced into the sub-iteration. It is well known that the Aitken's accelerator works properly if the target numerical sequence has a convergence. In this study, the target field to be found the convergence is set as the force field. The corrected force field is obtained using 3 estimations and the Aitken's algorithm during the sub-iteration process. Thus: The corrected force field F tn correct and resulting displacement field are used as a converged solution at t = t n . The overall procedure of the strongly coupled CFD-FEA at time t = t n in this study can be summarized as follows (see also 1. Compute the CFD where changing the fluid mesh from X tn-Δt correct to X tn predict and then obtain F tn 1 . 2. Convert the force field F tn 1 to the external pressure field (P ex,tn 1 ) on the FE model via a mapping process. 3. Compute the FEA where changing the external pressure field from P ex,tn-Δt correct to P ex,tn 1 and then obtain the displacement field u tn 1 . 4. Transfer calculated displacement field u tn 1 to the CFD mesh and then update the fluid mesh (X tn 1 ) via a mesh morphing process.

k = 2 or 3 (the corrector):
5. Compute the CFD where changing the fluid mesh from X tn-Δt correct to X tn 1 or X tn 2 , and then obtain F tn 2 or F tn 3 . 6. Obtain the external pressure field P ex,tn 2 or P ex,tn correct via a mapping process. 7. Compute the FEA where changing the external load field from P ex,tn-Δt correct to P ex,tn 2 or P ex,tn correct , and then obtain the displacement field u tn 1 or u tn correct . 8. Transfer calculated displacement field u tn 1 or u tn correct to the CFD mesh and then update the fluid mesh (X tn 2 or X tn correct ) via a mesh morphing process. 9. (k = 2): Go back to step 5.
(k = 3): Compute the CFD where changing the fluid mesh from X tn-Δt correct to X tn correct , and then return to step 1 and advance the time (t = t n + Δt).
As for step 1, X tn predict should be somehow defined. In this study, X tn predict is decided based on the linear extrapolation of the deformation from the previous time step, i.e., X tn predict is simply determined by extrapolating X tn-2Δt correct and X tn-Δt correct . The decomposition of the deformation from the FEA is also applied as with the case of the weakly coupled method. The interpolation of the pressure fields between consecutive coupling time steps is conducted by using Eq. 5 as well as the weakly coupled method. The damping ratio applied in the FEA model is increased at the early stage of the coupling in the same manner. Coupling time step size Δt is held to 0.019 s in full scale as well as the weakly coupled method. When the strongly coupled method is adopted, the solution was successfully proceeded even when α = 0.3 is applied. This may be due to a fact that the artificial added mass effect is mitigated by means of the sub-iteration processes.

Validation of two-way coupled method
In this chapter, a series of validations of the present weakly and strongly coupled CFD-FEA results is provided. For comparison, the experimental measurements from Ref. [17] are referred. The CFD and FEA calculations along with the data transfer and mapping computations are conducted using a workstation with Intel Xeon E5-1680 processor and with 8 cores parallel processing. Elapsed wall clock time for computing each method per physical 10 s (in full scale) is summarized in Table 4. As the data transfer and mapping processes are included in the weakly coupled method, it takes more total computational time than the one-way coupled method. Much longer computational time is needed in computing the strongly coupled method, since the sub-iteration processes are further included.

Rigid body motion
The RBMs of the ship evaluated from the respective numerical methods are compared with the experimental result. Time series of the pitch and heave motion of the ship among the experiment, the rigid body CFD results, and the two-way coupled results are compared in Fig. 11. Results from the rigid body CFD are presented after 25 wave cycles to exclude the transient part. For two-way coupled methods, the rigid body motions are derived from CFD computation results at each coupling time step. The results of the initial stage of the twoway coupled calculations during which the damping ratio of the FE model is increased are excluded in the figure. In addition, coefficient of correlation values of the RBM from numerical methods against the experiment are summarized in Table 5. All the numerical results represent the high correlation to the experiment.
As regards the RBMs from the weakly coupled method, plotted by green solid lines, overall agreement with the experiment or the rigid body CFD can be found. Unstable fluctuations are found around 20-30 s in those peak values and periods, as seen from the figures, which may be due to a fact that transient parts remained in the presented results. However, it should be noted that a convergence study in terms of the coupling time step size Δt will be necessary, as the accuracy of the weakly coupled method may be influenced by the coupling time step size to a large extent [24].
The RBM time series evaluated from the strongly coupled method are also comparable to those from the rigid

Vertical bending moment
Next, the VBMs evaluated from the respective methods and the experiment are compared. The VBMs are decomposed into the normal wave-induced component and the whipping components by applying the band-pass filter (BPF) to the time series results. Cut-off frequency of the BPF to discriminate these components is set at 0.3 Hz, considering the wave period and natural frequency of 2-node vibration on each model. A comparison of the wave-induced VBMs at the mid-ship section (SS5.5, see Fig. 1) is shown in Fig. 12a. Values in the vertical axis are plotted in a non-dimensional manner. In the result, positive value of M v_wave indicates the hogging moment. As seen from the figure, overall agreement of the wave-induced VBMs derived from the weakly and strongly coupled methods are found to agree with the experimental one, as well as the one-way coupled result. In addition, coefficient of correlation values of the VBM from the numerical methods against the experiment are summarized in Table 6. All the wave-induced VBM from the numerical results represent the high correlation to the experiment. On the other hand, the coefficient of correlation values represent a large variation in the case of the whipping VBM, especially the significant low correlation is found for the weakly coupled method result.
A comparison of the whipping VBM time series at SS5.5 is shown in Fig. 12b. Since severe weather/operational conditions with the wave height 10 m and with the full service speed are adopted in this study, the peak amplitudes of the whipping are remarkably large. When the weakly coupled results are focused on, the weakly coupled method underestimates the whipping components than the experiment. This fact may imply that the slamming impact force is not properly simulated. Since the CFD calculation in the present weakly coupled method has used the under-relaxation factor α = 0.05 in SIMPLE scheme, which may be too small according to conventional cases [18], thus the pressure temporal values might not be calculated appropriately. Meanwhile, the instability issue of the weakly coupled solution is found for the under-relaxation factor α = 0.3 as mentioned in Sect. 3.2.1. Decreasing the coupling time step Δt is considered to be one of the efficient measures [24] when the CSS-based procedure is used, and then, a convergence study in terms of Δt with α = 0.3 will be necessary. Nonetheless, it would require quite a few computational efforts to find decent coupling time step for the weakly coupled method, and thereby a benefit of the weakly coupled method would be reduced. Further careful investigations to solve such dilemma will be necessary if the weakly coupled method is used for the whipping evaluation.
When the strongly coupled result is focused, see Fig. 12b, the whipping component agrees well in terms of the peak amplitudes of the moments with the one-way coupled method or the experimental results. Regarding the whipping amplitudes, those from the strongly coupled method are comparable to the one-way coupled results. The whipping components in the one-way coupled results oscillate about 0.87 Hz, which is close to the natural frequency  of the 2-node vibration of the FE model in dry condition (0.894 Hz). Meanwhile, the frequency of the vibration of the whipping components from the strongly coupled results is observed around at 0.58 Hz. The reduction of the vibration frequency is owing to the added mass effect with regards to the elastic deformation. In fact, the natural frequency of the 2-node vibration of the experimental model in the wet condition was 0.699 Hz, while it was 0.946 Hz in the dry condition, which is comparable to the present simulation result. From these results, one may conclude that the present strongly coupled method can give a quantitative estimation on the whipping VBM considering the added mass effect. The wave elevation near the ship hull, point 'P' in Fig. 13 (7 m away from hull surface), is obtained from the present strongly coupled method, and then, the frequency is checked. The wave frequency from the strongly coupled method is compared with the rigid body CFD in Fig. 14. One may find the increase of wave amplitude around 0.3-0.6 Hz to some extent, which may be attributed to the radiation wave induced by elastic vibration, but its magnitude is still small. Thus, the wave damping may not influence the result so much in the present case. Note that as the present study is limited to only one wave condition, further case study will be necessary to consistently evaluate the wave damping effect.

Water pressure
Finally, local water pressure values on the hull are focused. Time series of the pressure at PS1 (bow flare) and PS2 (bottom at SS4.5) in Fig. 1 are compared in Fig. 15, among the experimental measurement, the rigid body CFD which is used for the one-way coupled method, and the strongly coupled results. In these figures, hydrostatic components are excluded. From Fig. 15a, both the numerical results have well captured the water impact pressure peak comparing with the experiment. In the experimental result, the second pressure peaks follow the first impact pressure, at t = 12 s and t = 21 s. These peaks are considered to have been caused by the 2-node elastic vibration of the model. One can find from Fig. 15a that such vibratory components are not reproduced using the rigid body CFD calculation. For the strongly coupled results, on the other hand, the vibratory components have been successfully reproduced. Next, the results on the pressure at the bottom point PS2, Fig. 15b, are discussed. One can find from the experimental result that the pressure history includes prominent vibratory component about 0. 6-0.8 Hz. This oscillation is attributed to the 2-node elastic vibration of the model, since the natural frequency of 2-node vibration in wet condition was measured at 0.699 Hz. This vibratory component, hereinafter referred to as the hydroelastic component, explains the discrepancy between the rigid body CFD and the experiment. When the strongly coupled result is looked at, on the other hand, the amplitude of the hydroelastic component agrees with the experimental result. The frequency of the hydroelastic component, about 0.5-0.6 Hz, in the strongly coupled result, is a little smaller than the experimental result, as it may be attributed to the difference of the structural model between the present FE model and experimental model. Very high-frequency vibrations are occasionally observed in the strongly coupled result, e.g., at around t = 16 s in Fig. 15b. As the convergences on the force field have been only attained as per Eq. 6 in the sub-iteration stages, it is estimated that rigorous convergence between force and displacement fields might not have been attained at some coupling time steps. To prevent this, it is recommended to increase the number of sub-iterations, which in turn leads to an increase of computational efforts, or to adopt other coupling techniques, e.g., the relaxation method [32]. Nonetheless, it can be concluded that the present strongly coupled method sufficiently captured the hydroelastic behavior of the water pressure.

Formulation of DBM and US
Using the present strongly coupled CFD-FEA results, the DBM and hull girder US are evaluated. In this study, an interaction formula between the VBM and DBM given by Amlashi  For simplicity, we assume that the hull girder US is governed by the longitudinal stress on the outer bottom panel [34], and then, the ultimate VBM and DBM are assumed to follow the beam theory. The longitudinal stress on the outer bottom panel when the hull girder reaches US is denoted as σ u . Given the section modulus of the target cross section Z v , M u,v can be expressed as follows: The sectional modulus of the ship at the mid-ship section is calculated to be Z v = 50.71 [m 3 ].
To estimate M db and M u,db , the formulae given by Tatsumi et al. [34] is employed. Let us assume that M u,db is considered to be the bending moment when the longitudinal stress on outer bottom panel reaches σ u under pure bending of a unitary double-bottom structure. According to Tatsumi et al. [34], when a unit cross section, section C see in Fig. 16, is taken from overall cross section of the double-bottom structure, M u,db can be calculated as follows: where I c and h c are the second moment of area and the height of the neutral axis from the outer surface on the section C, respectively. The sectional modulus of the section C is calculated to be Z c = 0.214 (m 3 ). The DBM M db is calculated as follows: where k p is a spring constant on the partial bulkhead [35], P c is the container load, q is the water pressure applied on the bottom, l db is the longitudinal length of the double-bottom between the partial bulkhead and watertight bulkhead, and k w is a spring constant on the watertight bulkhead. Note that q is derived from the present strongly coupled CFD-FEA results. For more detailed derivation of Eq. 10, Tatsumi et al. [34] is to be consulted.

Results and discussion
The hydrostatic, wave-induced, and hydroelastic components in the VBM and DBM at SS4.5 are evaluated from the strongly coupled CFD-FEA result and are compared in Fig. 17. The time span for the evaluation is 10-30 s. All the results in the figures are non-dimensionalized by M u,v and M u,db , respectively. σ u = 280 (MPa) is adopted to derive M u,v and M u,db . Regarding the VBM, one may find that the whipping component represents a large contribution to the M u,v around 20 s, up to 53% of M u,v . As to the DBM, the hydrostatic component M db_s accounts for the majority of the total DBM, i.e., 9% of M u,db . On the other hand, it should be noted that the proportion of the hydroelastic component M db_hyel has appeared at non-negligible level, which is almost at the same level with the wave-induced component M db_wave .
According to Eq. 7, L v , L db , and L comb are calculated and compared in Fig. 18 Fig. 19. It is found that the DBM component reduces the ultimate strength by about 11%. Similar tendency on the DBM contribution to the ultimate strength was reported in a precursory research [2], which observed the 18.3% reduction of the ultimate strength due to the DBM thorough the nonlinear FE analyses. Note that this reduction rate may change depending on the structural design of ship, assumed ship operation, etc. Finally, magnitudes of the hydrostatic, wave-induced, and hydroelastic components in M db at t = 25.43 s are broken down in Fig. 20. The largest component is the hydrostatic component representing 80% of the total DBM, as the full loading condition is assumed in this study. As regards the wave-induced and hydroelastic components, respective proportions are 11% and 7% of the total DBM. It can be interpreted that the hydroelastic vibration in the water pressure contributes the total DBM, eventually the US, to some extent. Attention should be paid in evaluating the DBM other than full loading condition, as the hydrostatic DBM may be decreased and thereby the proportion of the hydroelastic component may increase instead. In fact, it has been pointed out that the hydrostatic loads in ships have a significant statistical variation [36,37]. Since the present hydroelastic component in the DBM has been evaluated only by the water pressure fluctuation due to the whipping, an effect from the inertia force due to the whipping induced vertical displacement is not taken into account. According to the results presented in Ref. [16], this effect might also magnify the hydroelastic component in the DBM in the same phase. The two-way coupled CFD-FEA considering the flexibility of the double bottom will be needed to assess this effect quantitatively. Further research effort will be necessary to reach a general picture of the effect of the hydroelastic vibration on the DBM and subsequently the hull girder US.

Conclusions
In this paper, to establish a consistent evaluation method for the combined VBM and DBM, two-way coupled methods between the CFD and FEA are developed and validated. A 3D container ship model is used for evaluating the VBM and DBM based on the methods. Then, an investigation is made into the effect of the hydroelastic component of the DBM on the hull girder US. The followings may be concluded: 1. Good accuracies in terms of the RBM and the waveinduced components in the VBM and local pressure are confirmed when the present weakly coupled method is adopted. Discrepancies in the whipping vibrations arise from the poor representation of the added mass with regards to the elastic deformation. Further convergence studies in terms of the coupling time step size in the weakly coupled method will be necessary while using proper under-relaxation factor in SIMPLE scheme. 2. The present strongly coupled method gives better agreement with the experimental results in terms of the RBM, VBM, and local pressure, including the whipping vibration with the added mass with regards to the elastic deformation. The vibratory or hydroelastic component observed in the water pressure at the ship bottom also compares well with the experimental result. 3. The DBM evaluated from the present strongly coupled method reduced the hull girder US by about 11% in the considered case. Note that this value may change depending on the structural design of ship, assumed ship operation, or the location to evaluate the DBM [38]. 4. The hydroelastic component accounts for 7% of the total DBM. This indicates that the hydroelastic component in the DBM is non-negligible for accurate evaluation of the reduction of the hull girder US due to the DBM. The present strong coupling of the CFD and FEA is considered as a promising method for a consistent estimation of the hydroelastic effect on the VBM and DBM in a ship. Further investigation seems necessary in the future to clarify the effect of the DBM on the hull girder US taking the hydroelastic effect into account. Moreover, the extreme value prediction of the combined VBM and DBM with consideration of the hydroelastic vibration should be made in the future, as a follow-up work of our recent study [15].
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.