Efficient FORM-based extreme value prediction of nonlinear ship loads with an application of reduced-order model for coupled CFD and FEA

In this paper, a method for predicting the extreme value distribution of the vertical bending moment (VBM) in a ship under a given short-term sea state is presented. To predict the extreme value distribution of the VBM, the first-order reliability method (FORM), by which the most probable wave episodes (MPWEs) leading to given VBMs are identified, is introduced. The coupled computational fluid dynamics (CFD) and finite-element analysis (FEA) are used to provide the high-fidelity numerical solutions for the wave-induced and whipping components of the VBM. Then, a Reduced-Order Model (ROM) which can yield the predictions equivalent to the coupled CFD–FEA results in a relatively short time is developed. The ROM is incorporated into FORM to identify the MPWEs, in lieu of the coupled CFD–FEA. The accuracy of the ROM is verified by comparing with the coupled CFD–FEA results under identified MPWEs, in terms of both the wave-induced and whipping VBM. Then, a series of tank tests using a scaled container ship is conducted. In the first series of the tests, the VBM measurements under the MPWEs identified from the FORM-based approach using the ROM are made, to validate the accuracy of the ROM. The extreme value distribution of the combined wave-induced and whipping VBM is also measured by performing the second series of the tests, in random waves. The validity of the FORM-based extreme value prediction using the ROM is investigated by comparing with the second series of the tests.


Introduction
For ensuring structural safety of ships, a consistent estimation of the extreme value distribution of structural response is of particular importance. The long-term prediction method for any structural responses has been widely applied to determine the extreme response at various probability levels [1,2]. One can identify the most critical short-term sea state by referencing the contribution factor of each shortterm sea state to the long-term predictions. Such long-term prediction approach was extended to evaluating hydroelastic response of container ships, in particular, the slamming induced whipping response [3] or fatigue damage accumulation [4]. Although uncertainties may exist in the assumptions made for the structural response predictions [5], e.g., uncertainties due to operational speed, the effectiveness of determining the critical short-term sea state has been widely accepted.
Given the critical short-term sea state, the next challenge is to predict the extreme structural response in the short-term sea state including the various nonlinearities as precisely as possible. If one can ignore the computational costs, nonlinear time-domain simulations such as computational fluid dynamic (CFD) can be performed under the random irregular waves. However, it is impractical due to large computational burdens if the simulation method adopted is of high fidelity. Hence, the classification societies have attempted to provide the 'Design Regular Wave' concept, so that the extreme response level and the design load are approximated in a straightforward manner [6]. However, the accuracy of this approach is open to doubt due to the insufficient rationality in assuming the Design Regular Wave. To provide more accurate prediction, Adegeest et al. [7] and Fukasawa [8] suggested the estimation methods for predicting maximum value of structural response under irregular wave condition. Their methods are known as the most likely extreme response (MLER) method and the design irregular wave (DIW) method, respectively, from which conditional irregular wave trains are derived based on the linear superposition of response functions of structural response of interest under component regular waves.
From a view point of risk management of ships, on the other hand, it is highly necessary to evaluate the failure probabilities under the sea states for a given limit state. An effective method for predicting the probability of exceedance (PoE) was given by Der Kiureghian [9] based on the first-order reliability method (FORM). One can evaluate the extremes linked with the structural reliability or the PoE using the FORM. This FORM-based approach also enables us to identify the most probable wave episodes (MPWEs) associated with the maximum responses linked with a failure probability for a given linear/nonlinear failure function. An effectiveness of applying the FORM-based approach to extreme wave-induced vertical bending moment (VBM) prediction was demonstrated by Jensen [10]. Further extension of this work was conducted for evaluating VBM with nonlinearity taken into consideration, i.e., whipping effect of a flexible ship, by leveraging the nonlinear strip theory [11]. An application of the FORM-based approach to extreme value problems of structural responses subjected to combined loads was performed by Iijima et al. [12]. From these studies, it was indicated that the FORMbased approach can provide a strong tool for determining the MPWEs even when the structural response is nonlinear or is subjected to correlated load patterns.
Aiming at a derivation of the MPWEs for the extreme VBM by FORM, an accurate estimation of hydrodynamic forces including the slamming impact loads and subsequent whipping vibration is necessary. Given the noticeable advances of the computational performance in recent years, the CFD is considered to be one of the promising methods for evaluating them. Up until now, several research efforts have been spent for evaluating hydrodynamic forces under head sea conditions with the CFDs [13,14]. The expansion of the applicability of the CFD towards the whipping evaluations was carried out by integrating the structural model into hydrodynamic calculations [15,16]. Our research group has also employed the CFD to derive hydrodynamic forces on ship along with slamming impact loads and the whipping response. We developed a coupled method making use of the CFD and the finite-element analysis (FEA) in our preceding studies [17,18]. It was demonstrated that the coupled CFD-FEA approach can give a consistent estimation of the slamming impact force, whipping response, and the local bending response [19].
Several attempts to apply the CFD to the extreme VBM predictions have been made. Oberhagemann et al. [20] proposed a method based on the Monte Carlo simulations (MCS) and the CFD for small response levels while applying the Conditional Random Response Waves (CRRW) method [21] for large response levels. However, their method necessitated huge computational efforts in case of conducting long-term extreme value predictions. From a viewpoint of practical use, more efficient method to minimize the computational efforts while keeping the accuracy is needed. Seng et al. [22] adopted the open source CFD code OpenFOAM and the Timoshenko beam model to the extreme VBM prediction by FORM. They utilized the model-correction-factor (MCF) approach [23] in which the critical wave episode associated with the extreme response level is predicted using the nonlinear strip method (predictor step); then, the extreme response is corrected using the CFD-based method under the predicted wave episode (corrector step). A noteworthy point of this work is that they overcame an issue concerning the computational efforts in predicting extremes based on a combination of the CFD-based method and FORM, by employing the MCF approach. However, their study also posed a problem in evaluating the VBM with whipping component, which was fundamentally attributed to the discrepancies on physical model for slamming impact and whipping response assumed in the predictor and corrector methods. Our research group has also tackled with the extreme value prediction problem based on the coupled CFD-FEA method using the MCF approach [24,25]. It was reported that although the MCF approach is valid in predicting extreme values of the wave-induced VBM [24] or combined wave-induced VBM and local bending moment [25], a drawback of the approach was highlighted in predicting the extreme VBM with whipping component, which is consistent with the report by Seng et al. [22].
The above-mentioned background motivated our research activities regarding development of prediction method for extreme whipping responses based on the coupled CFD-FEA. To this end, an efficient alternative to the previous approaches, by which accurate prediction of the whipping response from the coupled CFD-FEA can be attained, is vital. During the last few decades, the reducedorder models (ROMs) as a surrogate for a more high-fidelity model have been proposed in several research fields [26], to represent high-fidelity numerical methods. Research efforts for devising ROMs were spent towards problems on fluid dynamics [27], structural responses [28], or fluid-structure interactions [29], for instance. Although ROM-based approach covers a wide variety of research areas, the fundamental concept of ROM development is common; substitute the high-fidelity models with a low-dimensional physic-based model, such that the dominant behaviors of interest can be reproduced via near real-time computations. A combination of the ROM and FORM (ROM+FORM) may be a potential alternative to the existing approaches for the extreme value predictions if the wave-induced and whipping VBMs are realized using simple mathematical formulations.
In this paper, a new approach for extreme VBM prediction considering whipping effect of a container ship, which is equivalent to a combination of the coupled CFD-FEA and FORM (CFD-FEA+FORM) is proposed. A pre-calibrated ROM so as to reproduce the coupled CFD-FEA results in an inexpensive manner is used for the extreme value predictions by FORM. Two distinct ROMs are proposed for the wave-induced and whipping components of VBM. Response predictions in terms of the wave-induced VBM are made utilizing the transfer function (TF) of the wave-induced VBM as the ROM. Rough screening of the TFs of the waveinduced VBM is carried out by leveraging the nonlinear strip theory, and then, the TFs are corrected based on the coupled CFD-FEA calculations. As an ROM for the whipping VBM, a simplified model based on a single-degree-of-freedom vibration system is employed. A classical momentum theory is applied for approximating the slamming impact force, and then, the whipping VBM is calculated by leveraging an impulse response function, along with some corrections for the coupled CFD-FEA. The present ROM+FORM approach is verified by comparing with the coupled CFD-FEA results under identified MPWEs. In addition, a towing tank experiment using a scaled model of a container ship in head sea conditions is newly conducted for a series of validations. A validation study of the present ROM is then carried out by comparing with the experiment under the MPWEs. Following this, the PoEs corresponding to extreme values of the wave-induced and whipping VBM are measured from the experiment, and these are compared with the prediction results based on ROM+FORM. Consequently, an insight on the appropriate evaluation of the stochastic characteristic of the whipping VBM is highlighted.

General
A fundamental theory for predicting the extreme values and MPWEs proposed by Jensen [11] is adopted in this study. When one considers linear, long-crested irregular waves, the free-surface elevation can be represented as a superposition of N discrete harmonic wave components: where S(ω i ) denotes the wave spectrum, ω i is discrete frequencies, and dω i is the increment between discrete frequencies. u i and ū i are the independent and Gaussian distributed stochastic variables. In this study, the ISSC wave spectra within ω i range of 0.3-1.5 [rad/s] are adopted. The number of discrete harmonic wave components N is taken to be 100. The frequency interval dω i is set to be non-equidistant to avoid repetition in waves. Let any structural response at the target time t 0 under the irregular wave train given by Eq. 1 be r t t 0 |u 1 , u 2 , ⋯ u N ,ū 1 ,ū 2 , ⋯ū N . In this study, response r t denotes the VBM estimated using any time-domain analyses, e.g., nonlinear strip theory, the coupled CFD-FEA, etc. According to FORM, the design point ( u * i ,ū * i ) is defined as a point on the limit surface, with the shortest distance to the origin, see Fig. 1. For the target extreme VBM M R , the limit state function g may be given as The reliability index β is obtained as a constrained optimization problem: The MPWE leading to the given M R is eventually predicted by assigning the design point ( u * i ,ū * i ) to Eq. 1. Once the reliability index β is given by Eq. 3, the PoE of M R may be given analytically as in Eq. 4 [30]:

Design point detection
An iterative scheme is required to detect the design point based on FORM if some nonlinearities, e.g., whipping  effect, are relevant in the target structural response. The general procedure of the iterative scheme is described in Fig. 2. At the first step of the process, combinations of u i and ū i are generated according to the Box-Muller algorithm, assuming that the initial design point is zero. In the subsequent processes, combinations of u i and ū i are randomly generated around the tentative design point ( u * i ,ū * i ), which is detected from the previous step, followed by updating ( u * i ,ū * i ). The limit state function g is approximated as a linear function of u i and ū i around the tentative design point, using the polynomial based Response Surface Method (RSM): Note that to evaluate the unknown coefficients c i and c i and to estimate the response surface which approximates the true limit surface around the tentative design point ( u * i ,ū * i ), substantial number of estimations of the limit state function g under various combination of u i and ū i , at least 4 N + 2 estimations, is needed. To solve Eq. 3, a Lagrangian function f is introduced as below with the approximated limit surface g = 0 as a constraint: Equation 3 may be readily solved by applying Lagrange multipliers method. It reduces to a set of simultaneous linear equations. As shown in Fig. 2, such updating processes of the design point are repeated until the convergence of design point will be found (inner iteration).

Outer iteraƟon
Even if the design point is found from the inner iteration schemes, it should be noted that the VBM under the estimated design point does not necessarily indicate the extreme value at the target time t 0 , in case of the target VBM includes the nonlinearities, e.g., sag-hog asymmetry and whipping after slamming impact. To obtain the true design point, further iterative scheme (outer iteration) is implemented during the detection processes, see Fig. 2. In the outer iteration, the target time t 0 will be updated depending on the occurrence time of extreme VBM, M e . The outer iteration process is repeated until the coincidence of the occurrence time of M e (i.e., t 0 ′ = t 0 ) and the convergence of M e are found. The convergence criterion, ε in Fig. 2, is set to be 0.001 in this study.

Subject ship and experimental setup
In this study, a series of towing tank tests using a scaled model of a recent container ship is conducted. The experimental results are used for validating the accuracy of the ROM + FORM approach, in terms of VBM value itself and PoE levels of VBM.

Ship model
A segmented scaled model of a recent container ship is used for a series of towing tank tests. The experiment using this model is carried out in the 150 m towing tank of National Maritime Research Institute (NMRI). The principal particulars of the model are described in Table 1 comparing with the prototype full-scale ship. The whole ship hull is 2.838 m in length, 0.428 m in breadth and 0.240 m in depth, and is divided into two rigid bodies amidships. Its scale ratio is assumed to be 1/100, see Table 1. As a result of ballasting and leveling of the model during the setup process, there are discrepancies of the total weight and draft between assumed full-scale ship and the scaled model, i.e., 11.6% heavier than the full-scale ship. The numerical models explained later are modeled based on the total weight of the scaled model.
Fundamental strategies for designing the model may be found in our previous works [31][32][33]. The segmented model consists of two rigid bodies each being connected with a hinge type device which also gives the flexibility to the structure. The bending rigidity is adjustable by changing the cross section of specimens equipped at the connecting section. This time, a beam shape specimen has been designed and installed amidships (see Fig. 3). The specimen bends and bears shear force and the hull girder deforms around the hinge when the hull girder is subjected to the VBM. The measurement of the VBM is conducted from the axial loads (L 1 , L 2, and L 3 ) measured via load cells (LUX-B-ID, Kyowa Electronic Instruments Co., Ltd.) placed at the cross section (see lower figure in Fig. 3). Given the measured axial loads and the vertical height of load cells H l , the VBM in Sect. 1 is calculated by [33] Here, the height H l is 0.130 m. The cross section of the specimen is called "trough type" [33] which aims to reproduce the realistic relationship between the bending moment and rotational angle of a ship hull girder, see Fig. 4. The scantling of the specimen is listed in Table 2. The natural frequency of 2-node vertical vibration mode of the model is measured in wet condition. Hammering tests to the model

Measurement
A series of towing tank tests using the ship model is conducted under irregular head seas and pre-determined MPWEs. Two short-term sea states described in Table 3 are selected for measurement. The forward speed of the model is zero knot and eight knots (in full scale) for State 1 and State 2, respectively. For State 1, the whipping effect is not expected to occur, whereas it is expected for State 2. Time histories of the irregular waves are prepared according to Eq. 1 by varying Gaussian distributed random numbers, u i and ū i . Each one irregular wave test lasts 5 min for State 1 and 3.5 min for State 2, respectively. It is repeated 80 times for different sets of u i and ū i to measure the extreme VBM statistics. These correspond to the total measurement periods of 66.7 h and 46.7 h at full scale, see Table 3. The PoEs of the VBM under State 1 and State 2 are evaluated from the measured VBMs. To count the individual peaks over the measurements, the zero-upcrossing periods are to be detected first from the wave-induced components (WF) of the VBMs. The waveinduced components from the measurement (WF) are derived by applying the band-pass filter (BPF) to the measured (with whipping, WF + HF) VBM. The cutoff frequency of BPF is set 3.0 Hz (0.3 Hz in full scale) to filter out the 2-node vibration components. Figure 5 provides an example. Once the zero-upcrossing periods are detected, the peaks of the individual waves are detected as the maximum (or minimum) value in the zero-upcrossing period. The number of counted peaks is 16538 and 13721 for State 1 and State 2, respectively. Figure 6 shows a snapshot of the installed model in the towing tank. The model is connected to the towing carriage via two towing rods. The towing test apparatus used in this study was not able to constraint the roll motion. Therefore, ship motions other than heave, pitch, and roll motions are constrained. The setup of the model under irregular waves on State 1 is schematically, as shown in

MPWE generation by wave maker
In the experiments, MPWEs are generated using a plunger type wave maker. To generate the target MPWEs, one should prescribe the vertical displacement of the wave maker which in turn the target wave fields can be successfully reproduced. According to Ref. [34], the relation between the wave elevation history η e and the wave maker motion can be given by where v denotes the velocity of wave maker, η I denotes the impulse response function of the wave maker, and θ denotes the angle of the wave maker, see Fig. 7a). Considering the two-dimensional water channel flow, η I is approximately given by where H e *, H I *, V*, and S* are the Fourier transform of η e , η I , v, and s, respectively. A series of calibrations of the MPWE generation by a wave maker is carried out using measured wave elevation via a wave height meter placed at X = 20[m]. The measured MPWE generated by the above-mentioned method is shown in Fig. 8, by a gray solid line. Measured wave elevation shows deviations from the target one in terms of the amplitude and phase. Since the impulse function η I given in Eq. 9 is an approximate expression based on the two-dimensional water channel assumption, some corrections are necessary by referencing to the wave-making signals according to the measured wave elevations. Given the measured wave elevation η m and its Fourier transform H m *, let the corrected vertical displacement of the wave maker s′(t) be (10) s where k(ω) is the weighting function in terms of ω, which is to be determined such that s′(t) induces a proper reproduction of target wave elevation. After several trials to determine the weighting function k(ω), a comparison between measured wave elevation and the target wave is shown in Fig. 8, by a black solid line. The amplitude and phase are successfully corrected via the calibrations. Thus, the MPWEs evaluated in this study are generated in the water tank using s′(t) in Eq. 12.

Coupled CFD-FEA
As a method for computing VBM using the high-fidelity numerical models, a coupling method between the CFD and FEA [19] is employed in this study. The CFD based on the finite-volume method (FVM) is used for the hydrodynamic force evaluations exerted on the ship hull. Fundamental methodologies of the CFD calculation conditions are compliant with our previous studies [19,24]. The commercial solver STARCCM + 13.06.012 [35] is adopted for the CFD computations. The CFD computations are conducted in the model scale, where the assumed ship length L pp is 3 m, and then, the hydrodynamic pressures over the discretized hull surface meshes are calculated. Hereafter, all the numerical results will be presented at full scale instead of the model scale. Figure 9 provides the solution domain of the CFD. By taking account of the symmetry of the problem, only half of the ship model cutting off at Y = 0 surface is used in the solution domain, while a symmetry boundary condition is set at Y = 0 section. The Euler overlay method (EOM) [36,37] is applied to the present CFD to reduce effects from the reflection of surface waves from the boundaries or the freesurface disturbance due to the ship motion. The width of the solution domain is set 0.8 L pp which may be enough to avoid reflected surface waves from the Y-direction disturbing flows near the ship [38]. The region of the damping zone in applying EOM is defined as is the case in Ref. [19].
Three CFD models are used to check the variability of numerical results. Commercial FEA solver in LS-DYNA ver_971R9.0 is adopted. A description of each CFD model is shown in Table 4. Amongst model 1 and model 2, the mesh resolution of the CFD over the free-surface region, Δx (horizontal direction) and Δz (vertical direction), has been changed. Model 3 has the same mesh resolution with Model 1, but the time increment during the CFD computation is decreased to 0.01 s. Elapsed wall clock time for computing the CFD per physical 100 s under 8 cores parallel computation are also described in Table 4.
The VBMs are calculated from the FEA computations. Since the VBM can be realized as a response of a beam subjected to the vertical forces, a 1D beam FE model is introduced for evaluating the hull flexibility. In this study, 20 beam elements of which overall length are equal to L pp are used. The mass distribution along the beam is taken from Ref. [19] based on that of the full-scale ship, and then, each mass of the elements is equally multiplied, so that the total mass of the model becomes equal to that of the experimental model. Young's modulus of the beam elements is adjusted, so that the natural frequency of the 2-node vertical bending mode of the model (see Fig. 10) becomes equal to that of the experimental model. In performing the FEA, 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 [18]. In coupling the CFD and the FE model, the one-way coupling between CFD and FEA [19] is adopted. The CFD mesh on the hull surface is discretized into 21 groups then the loads integrated over  the respective groups are exerted to the respective structural nodes. There may be small errors in terms of the interpolation converting from the loads on the three dimensional hull surface to the one dimensional forces. In this study, correction factor c l is introduced to balance the total force on the numerical models between CFD and FEA. It means: where f ZFEA,i denotes the vertical force acting on the ith FE node, f Z,CFD denotes the vertical force acting on the hull at coordinate x, and L i denotes the length of integral range. The vertical load distribution on a calm water surface is preliminarily calculated by the CFD, and then, the value of c l is determined, so that the sum of the vertical loads is balanced between CFD and FEA. In this study, the VBM amidships are targeted to be evaluated. The wave-induced VBM and whipping VBM components are discriminated via the BPF.

Reduced-order model for VBM
Here, a new ROM which gives the peak VBM equivalent to the CFD-FEA-coupled method with less computational effort is explained. The ROM is used for a series of extreme value predictions by FORM (ROM+FORM) aiming at the prediction of the design points corresponding to the CFD-FEA+FORM.

ROM for wave-induced VBM
A prediction of the wave-induced component of the VBM is made based on the TF and some corrections of peak values to account for the sag-hog nonlinearity of the wave-induced VBM. It should be noted that the solution time must be taken long enough for the simulation to reach steady states to evaluate the TF correctly. However, it would need quite a few computational efforts to obtain time-domain results from the coupled CFD-FEA under irregular wave. Hence, in this study, the nonlinear strip theory implemented in in-house code NMRIW-II [39] is used for obtaining the TF in lieu of direct computation of the coupled CFD-FEA. Let the wave-induced VBM calculated from NMRIW-II under the irregular wave expressed by Eq. 1, M wave,strip , be as follows: where R(ω i ) and T(ω i ) are the amplitude and phase of the TF, respectively. c strip means the correction factor for nonlinearity. Supposing that the nonlinearity of wave-induced VBM is relevant to the amplitude of wave elevation, let c strip be approximated by an nth degree polynomial in terms of peak wave amplitude η p : In this study, n = 6 is adopted. Note that although η p [m] is a dimensional value, it is treated as dimensionless temporarily in Eq. 17. The procedure to determine the TF and c strip is summarized as follows: 1. Generate random variables u i and ū i then compute NMRIW-II to obtain the wave-induced VBM. MPWEs, then a k in Eq. 16 is determined using the least-square approach, so that the extreme wave-induced VBM at target time t 0 from NMRIW-II is expressed by Eq. 15.
As for step 1, the wave elevation histories of 7500 s (for State 1) and 6500 s (for State 2) in full scale are used to obtain sufficient amount of data for the FFT. The amplitude spectrum of the wave-induced VBM from NMRIW-II under State 1 is shown in Fig. 11. Derived TF of the VBM, R(ω i ) and T(ω i ), calculated from NMRIW-II under State 1 is shown in Fig. 12. From Fig. 12, significant high local peaks of R(ω i ) are found in particular when ω is larger than 1.0 [rad/s]. These are attributed to the fact that the contribution of higher wave frequency components to the wave and wave-induced VBM amplitudes is small, as found from Fig. 11, but the round-off error by the FFT results in such a R(ω i ) with many local peaks. The local peaks should affect the quality of the MPWE if the local peaks are considered when obtaining the MPWE. To avoid this, the local peaks in the curve are smoothed out. In addition, since the high-frequency component of the wave-induced VBM is negligibly small, as found from Fig. 11, high-frequency (15) M wave,strip (t) After R(ω i ), T(ω i ), and c strip are determined, a further correction to M wave,strip is made to match with the coupled CFD-FEA results. According to our preceding study [24], the extreme wave-induced VBM from the coupled CFD-FEA may be predicted via a simple correction of peak amplitudes of the wave-induced VBM from NMRIW-II. Then, the wave-induced VBM is corrected by multiplying a factor c cfd , which is to be determined based on the coupled CFD-FEA results: where c cfd is defined as a function of η p , and M wave,CFD is the wave-induced VBM corresponding to the coupled CFD-FEA results. To estimate c cfd , several MPWEs are prepared by combining updated Eq. 15 and FORM, then c cfd in terms of η p is evaluated based on the coupled CFD-FEA results under the MPWEs. As a consequence, M wave,CFD is used as the ROM.

ROM for whipping component of VBM
Next, an ROM for whipping component of VBM is formulated. The first step to this end is to estimate the TF of ship motion from the CFD, viz. heave and pitch motions. These TFs can be evaluated in the same way with estimating R(ω i ) and T(ω i ) in Eq. 15, without considering the nonlinearity. Thus where h and θ are heave and pitch motion, respectively, and (R h , T h ) and (R θ , T θ ) are the TFs of heave and pitch motion, respectively. Heave motion h takes positive values along the Z-axis upward and positive values of pitch motion θ follow the rotation about Y-axis based on the right-hand rule, see Fig. 13. Given the heave and pitch motion histories, the time history of relative distance between the bow and the wave surface, η rel , is expressed as follows: where l bow denotes the distance between the center of gravity (CoG) of the ship and target bow section. Note that there is an upper limit to the magnitude of η rel , i.e., η rel is less than the deck height of the ship. η bow denotes the wave elevation from the calm water surface at target bow section, see Fig. 13, which can be calculated according to Eq. 21: where k i is the wavenumber, x bow is location of target bow section from amidships. The slamming impact force acting on the bow is approximated by the two-dimensional water impact assumption. According to the Karman's momentum theory [40], let the slamming impact force F imp on a two-dimensional wedge profile be as below: where the dot above the variable denotes a differentiation with respect to time t. M a denotes the virtual added mass, which can be derived from the following equation without considering the water pile-up effect [41]: where ρ is the fluid density, η rel -d denotes the height of the wedge beneath the water surface, and δ is the deadrise angle of the wedge against calm water surface, see Fig. 14. One can adjust the time of onset of slamming impact and its magnitude by arbitrarily changing d and δ.
Supposing that the whipping vibration can be realized by the oscillatory system with one-DoF. In this sense, the whipping VBM M whip is expressed as follows: where I is the impulse response function. c imp and τ imp are the correction factors for the magnitude and phase of slamming impact force, respectively. As the main contribution to the whipping VBM is 2-node vibration of the ship beam, which may be realized by the damped free vibration, I can be approximated using an analogy to one-DoF system subjected to an impact: where ζ and ω 0 are the damping ratio and the natural angular frequency of 2-node vibration, respectively. In this study, I(0) = 0 and ̇I (0) = −1 are adopted to express the impulse response function.
Once the ROM for the wave-induced VBM is determined according to the procedures mentioned in Sect. 5.1, an MPWE for arbitrary extreme wave-induced VBM can be estimated by combining with FORM. The TFs for heave and pitch motions are then determined from the CFD result under an MPWE. Following the coupled CFD-FEA computation, the whipping component of VBM is extracted by removing the wave-induced component using the BPF. Through a screening process of optimal values of d, δ, c imp , and τ imp , the ROM for whipping VBM is eventually completed. In this paper, d, δ, c imp , and τ imp are determined in a heuristic manner.

Results and discussion
The variability of the present coupled CFD-FEA method is first examined by comparing with the numerical results using different CFD models. Then, the accuracy of the present ROM + FORM approach on the wave-induced and whipping VBMs is verified by comparing with the coupled CFD-FEA results under the MPWEs. Finally, a series of validations is carried out by comparing the present ROM + FORM results with the experimental measurements.

Variability of Coupled CFD-FEA Method
Assuming that c cfd in Eq. 17 is equal to 1.0 on the ROM, i.e., M wave,CFD is assumed to be equal to M wave,strip , an MPWE which takes extreme wave-induced VBM M R of 3000MNm over State 1 at target time 600 s, named MPWE_S1a, is identified from the ROM + FORM computation. The wave elevations' amidships under the MPWE_S1a which are reproduced with the CFD models are shown in Fig. 15. In the figure, the black solid line denotes the target wave elevation given by Eq. 1. Although the Model 1 and Model 2 results represent slightly higher wave amplitudes on the troughs and crests around 600 s, all the present CFD models give the wave time histories close to the target wave elevation. When the results from Model 1 and Model 2 are compared, a quite small difference in the wave elevations is found despite the different mesh resolutions, see Table 4. This may indicate that further mesh refinement over Model 1 is not necessary for the solution accuracy. The result from Model 3 is close to that from Model 1 or Model 2, but the accuracy of the representation on the crest around 600 s has been slightly improved. Comparisons of the wave-induced VBMs among results from the coupled CFD-FEA results are shown in Fig. 16. A result from the ROM (c cfd = 1.0) is also plotted as a reference. From Fig. 16, one may find that all the CFD-FEA results resemble each other.
Coefficient of correlation values against the Model 1-based results in terms of the wave elevations and waveinduced VBMs are summarized in Table 5. Coefficient of correlation values is derived between 560 s and 620 s. It can be found that the variability of the wave-induced VBMs calculated by each CFD-FEA model is extremely small. Hereafter, Model 1 is used for the coupled CFD-FEA computations to calibrate and verify the ROM.

Wave-induced VBM
The correction factor for expressing the CFD-FEA results (c cfd , in Eq. 17) in the ROM is examined. Four CFD-FEA calculations are conducted for each short-term sea state to estimate the correction factor. Extreme hogging VBMs are targeted to determine c cfd . The estimated c cfd under State 1 and State 2 are plotted in Fig. 17 in η p terms. The magnitudes of c cfd are approximately 0.7 under State 1 and State 2, in the case of high-wave amplitudes. These results imply that the difference of the extreme VBM values between the coupled CFD-FEA and NMRIW-II is approximately 1.5 times. One possible cause of this discrepancy is that as the nonlinear strip theory adopted in the present, NMRIW-II calculations do not account for so-called memory effect [42], its accuracy of calculation might be poor in irregular waves. Using ROM incorporating c cfd and FORM, predictions of MPWEs associated with an extreme wave-induced VBM are  The present ROM+FORM approach is somewhat different from the known MCF approach [23], since the corrector stages using the coupled CFD-FEA results at obtained design points have not been included. To reach more rigorous design point corresponding to a combination of the coupled CFD-FEA and FORM, adopting the MCF approach, where the ROM+FORM and coupled CFD-FEA are used as the predictor and corrector, respectively, would be effective. Nonetheless, the present ROM + FORM was found to be able to predict already the design point of extreme hogging wave-induced VBM in a sufficient manner without the corrector stage, through pre-calibrations of the ROM with several coupled CFD-FEA results.

Combined wave-induced and whipping VBM
Next, the verification of the ROM + FORM approach for combined wave-induced and whipping VBM is made. State 2 is the subject sea state, as the whipping effect is negligibly small under State 1. The MPWE corresponding to Fig. 18b) is exploited to inspect appropriate values of d, δ, c imp , and τ imp in Eqs. (23) and (24). In this study, slamming impact forces are calculated at 6 cross sections on the bow (SS9.75, SS9.5, SS9.25, SS9.0, SS8.75, and SS8.5), and then, these sum is used for the whipping evaluation.  Fig. 19 along with the reliability index β. From the figure, it can be found that the reliability index β varies among each MPWE. A possible cause of such multiple design point problem is that several local minima/ maxima have been detected when the whipping component is included. Recently, the similar problem has also been reported in Ref. [43], where extreme parametric roll prediction was addressed by FORM. The authors of Ref. [43] dealt with this problem by classifying the design points as 'global minima' and 'local maxima close to global minima'. Such classification of the design points may be needed in the future, for evaluating appropriate reliability index β. The following verification of the ROM + FORM approach The present ROM could capture the peak VBM within 8% error in the presented cases; thus, it can be interpreted that the MPWEs equivalent to the CFD-FEA + FORM are sufficiently derived using ROM + FORM. Note that for reaching further rigorous design points, adopting the MCF approach with proper MCF would be ideal as well as the wave-induced VBM case. As the present ROM was found to be able to predict already comparable VBM to the coupled CFD-FEA both in terms of the amplitude and phase at the derived design points, it is expected that the MCF approach will also successfully work.

Comparison with towing tank test
This section is devoted to a validation of the proposed ROM + FORM approach by comparing with the towing tank test results. First, the accuracy of the VBM with the whipping vibration predicted by the present ROM is validated against the experimental result under the pre-determined MPWEs. The accuracy of the PoEs of the VBM with and without whipping is further validated by comparing with the experimental measurement under State 1 and State 2.

Validation under deterministic MPWEs
By applying the ROM + FORM approach, two MPWEs under State 2 are identified. The target extreme responses are the combined wave-induced (WF) and whipping (HF) VBM, while the target time is 150 s. Assumed M R are: 3900MNm (MPWE_S2a) and 4700MNm (MPWE_S2b). To examine the repeatability of the experiment, five trials are carried out on each case. Measured VBMs from these trials under MPWE_S2a and MPWE_S2b are shown in Fig. 21. A slight deviation in terms of the phase of VBM arisen from the whipping component, which might be due to the difference in the initiation time of the slamming impact, is found from Fig. 21. On the other hand, the variation in the amplitudes of the whipping vibration among each trial is found to be quite small.
Comparisons of VBMs between results from the present ROM and the experiment under MPWE_S2a and MPWE_ S2b are shown in Fig. 22. The fifth trial results are plotted as the experimental results. Dotted lines in the figures denote the results on the wave-induced component. From Fig. 22, one can find that the present ROM could well predict the wave-induced VBMs which are comparable to those from the experiment. As to the whipping component, the discrepancy of the amplitude of the whipping near target time (150 s) has appeared in the case of MPWE_S2a, see Fig. 22a). Regarding the result under MPWE_S2b, on the other hand, a good correlation between the experiment and the present ROM has been indicated in terms of the amplitude of the whipping, see Fig. 22b). As is mentioned in sub-sect. 6.2.2, it would be improved if the values of d, δ, c imp , and τ imp are determined depending on the target extreme VBM.

Validation of PoE by ROM + FORM
Finally, the PoEs of the VBM (with and without whipping) under State 1 and State 2 are compared between the experiment and the ROM + FORM approach. As already described in sub-sect. 6  the ROM + FORM computations in the case that nonlinearity becomes larger; particularly when the whipping component is included. The first approximation would be to adopt the minimum value of β and the associated design point. Therefore, in this study, 20 computations of ROM + FORM are conducted, and then, the minimum β are adopted to evaluate PoEs. Figure 23 provides a comparison of the PoEs under State 1 and State 2. The respective plots are linked with the value of β on the design points. The result for State 1 includes only the wave-induced VBM, while that for State 2 includes both the wave-induced VBM (WF) and the combined waveinduced and whipping (WF + HF) VBM. It is observed that the present ROM + FORM approach predicts the PoE of wave-induced VBM in an acceptable manner for both State 1 and State 2. The present ROM + FORM seems to have slightly underestimated the PoEs of the wave-induced VBM near 4000[MNm] under State 2. A correction to the present ROM may be necessary to represent the PoEs in the large wave-induced VBM range more accurately. Since the present ROM for the wave-induced VBM is constructed using the TFs based on the nonlinear strip theory, there is still room for improvement using the TFs based on the coupled CFD-FEA.
The discrepancy becomes larger when the VBM includes the whipping component (WF+ HF), in the case of State 2. Though the ROM + FORM results are presented assuming that d = 21.0[m], δ = 50.0[deg], c imp = 13.5, and τ imp = 0. 8[s] in the whipping ROM, there may be room for further optimization of these variables. As previously mentioned in sub-sect. 6.3.1, the screening of these variables may need to be made depending on the levels of M R . In addition, as is mentioned in sub-sect. 6.2.2, the classification of the multiple design points may be necessary to evaluate appropriate β. Further investigation about these points should be made for rigorous assessment of the PoE.
In the present ROM, only the bow flare slamming impact is accounted for. Thus, the stern slamming effect may need to be considered in the ROM with the same approach as the bow flare slamming. Moreover, the springing effect on the VBM might be non-negligible. As a reference, Fig. 24 shows an example of measured VBM time series from the experiment. From the WF + HF result in Fig. 24, high-frequency vibration components can be seen around 1300-1350 s despite the small wave height. The springing may have influenced to the PoE levels over the whole VBM region. One should consider the effect of springing together with whipping to capture more reliable PoEs. Further re-construction of the ROM addressing these causes may be pointed out as a future work.

Conclusions
In this paper, a new ROM is introduced and used with FORM for identifying the MPWEs leading to extreme hogging VBMs of a container ship, to approximate the extreme value predictions based on the CFD-FEA + FORM. The present ROM is formulated based on the hybrid use of the nonlinear strip method and the pre-calibrations using specific coupled CFD-FEA results. The ROM + FORM approach has been validated by comparing the tank test results. The followings are concluded.
1. The present ROM+FORM approach can consistently predict the design points for the extreme wave-induced hogging VBM equivalent to the CFD-FEA+FORM. 2. The present ROM+FORM approach can sufficiently predict the design points for the extreme hogging VBM equivalent to the CFD-FEA+FORM, even when the whipping component is included. Further application of the MCF approach will make it possible to reach further rigorous design points. 3. Through the comparisons between the ROM and the towing tank test results under pre-determined MPWEs, it turned out that the predicted VBMs from the present ROM are also comparable to those from the experiment, in terms of the wave-induced VBM and the amplitude of whipping VBM. It would be improved if the variables in the whipping ROM are determined depending on the extreme VBM levels. 4. The proposed ROM+FORM approach can consistently predict the PoEs of the wave-induced VBM. The results are comparable to the experimental results. There is still room for improvement in the large wave-induced VBM range using the TFs based on the coupled CFD-FEA. 5. The proposed ROM+FORM tends to underestimate the PoEs of the combined wave-induced and whipping VBM compared with the experiment. The underestimation of the PoEs may be partly attributed to the variable setting in the ROM or the multiple design points detected from the ROM+FORM approach. 6. Meanwhile, the present ROM needs to be extended to account for the high-frequency vibrations from the stern slamming and the springing effect. It can be done by following a similar procedure introduced in the present paper.