4400 TEU cargo ship dynamic analysis by Gaidai reliability method

Modern cargo vessel transport constitutes an important part of global economy; hence it is of paramount importance to develop novel, more efficient reliability meth-ods for cargo ships, especially if onboard recorded data is available. Classic reliability methods, dealing with timeseries, do not have the advantage of dealing efficiently with system high dimensionality and cross-correlation between different dimensions. This study validates novel structural reliability method suitable for multi-dimensional structural systems versus a well-established bivariate statistical method. An example of this reliability study was a chosen container ship subjected to large deck panel stresses during sailing. Risk of losing containers, due to extreme motions is the primary concern for ship cargo transport. Due to non-stationarity and complicated nonlinearities of both waves and ship motions, it is challenging to model such a phenomenon. In the case of extreme motions, the role of nonlinearities dramatically increases, activating effects of second and higher order. Moreover, laboratory tests may also be questioned. Therefore, data measured on actual ships during their voyages in harsh weather provides a unique insight into statistics of ship motions. This study aimed at benchmarking and validation of the state-of-the-art method, which enables extraction of the necessary information about the extreme system dynamics from onboard measured time histories. The method proposed in this study opens up broad possibilities of predicting simply, yet efficiently potential failure or structural damage risks for the nonlinear multi-dimensional cargo vessel dynamic systems as a whole. Note that advocated novel reliability method can be used for a wide range of complex engineering systems, thus not limited to cargo ship only.


Introduction
Extreme value prediction problems in engineering being frequent and challenging, especially when data is scarce (ISSC 2009;Wang and Wang 2021;Dulebenets 2022;Drummen et al. 2006;Abioye et al. 2021;Elmi et al. 2022;Chen et al. 2021;Jovanović et al. 2022;Parunov et al. 2015; European Maritime Safety Agency (EMSA) 2021).Hence, from a practical standpoint, it is essential to develop new, effective, and precise extrapolation approaches.Since container ships play a significant role in the contemporary marine transport business, this research is primarily concerned with their safety and reliability.The MSC Napoli, a post-Panamax container ship, experienced an engine room bulkhead failure in January 2007.An evaluation of the MSC Napoli's hull girder reliability at the time of the disaster was given by Parunov et al. (2015).The hull girder ultimate strength served as a representative of structural capacity, and the loads taken into account included still water bending moments, vertical wave bending moments, whipping bending moments.The load combination of whipping and wave bending moments was taken into account, 1st order reliability approach was used to calculate the likelihood of failure in Parunov et al. (2015).Second post-Panamax container ship, the MOL Comfort, also was broken down in June 2013 (Gafero Priapalla Rahim 2019; Falzarano et al. 2012;Su 2012;Ellermann 2008;BV, Bureau Veritas. 2012;Gaidai et al. 2023a).Both cargo ships were damaged as a result of overloading the main deck panel girder, giving them lower collapse strengths than other comparable vessels, even though they may not have been built and approved in accordance with best practice.Such significant events, especially those involving container ships, have rattled the maritime industry and demand a thorough investigation.Next, one can mention human error reliability aspect, within maritime transportation industrial sector (Bicen et al. 2021).Human error can occur for a variety of causes, including inadequate training and experience, poor communication, insufficient system monitoring, failure to learn from past mishaps, exhaustion, and stress.
The study's research object is a 245-m long 4400 TEU Panamax container ship (ISSC 2009).Motion sensors were installed on the 4400 TEU Panamax container ship during its trans-Atlantic journeys in 2010s years (Dulebenets 2022;Drummen et al. 2006;Abioye et al. 2021;Elmi et al. 2022).Captain made special measures to avoid the worst storms.4400 TEU dynamic ship panel stress sensor placement was carried out according to DNV container vessel standards and regulations (Chen et al. 2021;Jovanović et al. 2022;Parunov et al. 2015; European Maritime Safety Agency (EMSA) 2021; Gafero Priapalla Rahim 2019; Falzarano et al. 2012;Su 2012;Ellermann 2008).Onboard data sampling time step was dt = 0.025 sec, measured continuously over the whole voyage duration.For detailed ship parameters see (BV, Bureau Veritas 2012;Gaidai et al. 2023a;Bicen et al. 2021;DNV-RP-H103. 2011).Onboard measured data was not post-processed, namely raw dataset has been analyzed; all voyage-datasets were combined into one 1D (1-Dimensional) jointly representative time-series.Assessing extreme ship deck panel stresses and corresponding low probability levels is a challenging task, and only limited progress has been made in the past decade (Tomic et al. 2018;Sun et al. 2023;Gaidai et al. 2023bGaidai et al. , 2010a;;Andersen and Jensen 2014;Storhaug and Moe 2007).Numerical simulations being based on the underlying model assumptions, and being not able to capture all related wave and ship deck panel stresses nonlinearities.
When calculating the design loads under current ship regulations, ships are mostly treated as rigid.Cargo ships are actually flexible, and waves cause their hull structural vibrations.Critical loads and fatigue damage are exacerbated by these hull vibrations.Full scale measurements on Capesize bulk carriers provided evidence of major contributions to fatigue damage, but the impact on container ships overall reliability has not yet been thoroughly investigated.Large 294-m cargo ship sailing in North Atlantic has been used in this study to evaluate areal deck panel stress levels.Wave radar and wind data have been used to complement stress records amidships, and those measurements were evaluated aboard in real time by hull monitoring system to give captain valuable information.Warping, axial forces, horizontal and vertical bending moments-all contribute to increased stress levels and fatigue age.The difference between the overall fatigue damage, and the wave frequency damage is often used to assess fatigue damage, brought on by wave-induced vibrations.Wave-induced vibration effects depend on vessel size, eigen-periods, materials, presence of initial cracks.In order to decrease maintenance downtime throughout the course of vessel's service life, wave-induced vibrations effects should be taken into account during design verification for any large container vessel.Onboard measured data is vital to assess structural reliability, calibrate numerical tools, following design ship regulations.
Built in the year 2003, 4400TEU vessel followed contemporary design.DNV assessed onboard stress records from a 4000TEU container ship of a comparable size and age as part of AOG (Active Operator Guidance) research project.Full-scale measurements running in the North Pacific showed a contribution of between 39 and 46% depending on the season, in contrast to the 37% indicated by the head sea model experiments (Drummen et al. 2006).
Onboard 4400TEU vessel full scale measurements along with wave-induced vibrations analysis tools were kindly provided by DNV.
Figure 1 shows the 4400TEU cargo ship, used in trade between Germany, France, UK, Canada.Table 1 presents summary of the vessel's essential information.The vessel has the following notations: + 1A1, Container Carrier, EO, NAUTICUS (Newbuilding), ICE-C, W1-OC, DG-P, constructed to DNV class.The vessel has a typical cross-sectional shape but runs with a somewhat lower draft than other boats of this size.For instance, a lot of HT32 steel is utilized, with HT36 above the top deck.The hatch coaming plating is 65 mm thick, whereas the shear strake and top deck are built of 60 mm plates.The sensors SW and IS, shown in Fig. 5 were intended to be placed in the stiffener web's neutral axis as well as at the stiffener's momentary zero point for bending owing to lateral pressure.

Fig. 1 Example of loaded TEU container ship
In the next we will introduce novel reliability methodology in some details, following by application of Gaidai reliability method to the cargo ship onboard dataset.

Gaidai reliability method
Estimating structural system safety using traditional engineering reliability methods is often challenging (Wang and Wang 2021;Dulebenets 2022;Drummen et al. 2006;Abioye et al. 2021;Elmi et al. 2022;Chen et al. 2021;Jovanović et al. 2022;Parunov et al. 2015; European Maritime Safety Agency (EMSA) 2021; Gafero Priapalla Rahim 2019; Falzarano et al. 2012).The reliability of a complex structural system may be directly assessed by doing direct numerical Monte Carlo simulations, or by having sufficiently large measured dataset.Most of dynamic engineering systems are often too complex for both direct computational, and experimental methods.The authors have developed a novel measurement-cost-reduction approach for structural systems in an effort to reduce measurement costs, (Gaidai et al. 2023e, 2022k, 2022l, 2023f ).
Ocean waves are commonly thought to follow a stationary, homogenous ergodic random process throughout a 3-h storm.Think of a structure with several degrees of freedom that is subject to random stationary external loadings that are stationary in time, such as environmental wind and waves.The other option is to consider the process as   being reliant on a few outside elements, whose change over time may be referred to as an ergodic process.The dynamic system vector being composed of measured, or numerically simulated (X(t), Y (t), Z(t), . . . ) system key components X(t), Y (t), Z(t), . . .syn- chronously recorded over a sufficient (representative) -time lapse (0, T ) .The dynamic system 1D key component's global maxima, measured over the whole (representative) time lapse (0, T ) being denoted as . ., respectively.By large enough value of T authors mean large value in comparison to dynamic system relaxation and auto-correlation time scales (Gaidai et al. 2023a(Gaidai et al. , 2023b(Gaidai et al. , 2010a(Gaidai et al. , 2010b(Gaidai et al. , 2018(Gaidai et al. , 2016(Gaidai et al. , 2022a(Gaidai et al. , 2022b(Gaidai et al. , 2022c;;Sun et al. 2023;Andersen and Jensen 2014;Storhaug and Moe 2007;Balakrishna et al. 2022).Let X 1 , . . ., X N X be the temporally subsequent dynamic system key component's process within observational span (0, T ) and being either system's key loads, or key response on their own.Local maxima for other MDOF dynamic system components, such as Y (t), Z(t), . . ., being designated here as Y 1 , . . ., Y N Y ; Z 1 , . . ., Z N Z etc.It has been assumed that the system's key components local maxima being positive.Calculating the target failure/hazard probability/risk for dynamic systems is the current objective with being the target system non-exceedance (survival) probability for the critical values of the dynamic system's componentsη X ,η Y ,η Z ,…; ∪ denotes the logical unity operator «or»; and p X max T ,Y max T ,Z max T ,... representing joint PDF of the global system key component's maxima within(0, T ) .Due high dimensionality of most engineering dynamic systems, as well as limitations of the available dataset, it is not possible to explicitly estimate the later system's joint PDF, p X max T ,Y max T ,Z max T ,... in practice.In other words, dynamic system being considered having instantaneously failed (or entered in a state of hazard) at the time instants, when either X(t) exceedingη X , or Y (t) exceedingη Y , or Z(t) exceedsη Z , and so on, the fixed failure/hazard levelsη X ,η Y ,η Z ,… being unique for each single one-dimen- Occur- rence times of system key component's local maxima being denoted byt j .MDOF risk/ hazard/limit vector(η X , η Y , η Z , ...) , being introduced, consisting of hazard/failure limits of dynamic system's key componentsX, Y , Z, . . . .1D dynamic system's key response/load key component's local maxima being combined into a single merged synthetic (1) , according to the merged temporal system's vec- tort 1 ≤ • • • ≤ t N .In order to reduce limit/hazard values for each response/load compo- nent, the scaling parameter 0 < ≤ 1 being now introduced.MDOF limit/hazard vector η X , η Y , η z , ... can be now introduced, havingη The introduction of the unified system limit/hazard/risk vector η 1 , . . ., η N , with each component η j being eitherη X , η Y orη z , …, depending on which specific system key com- ponent supplied particular local maxima at the given time momentt j ; in case of simultaneous occurrence of several system key component's local maxima at the specific time instantt j , η j will be taken as a maximum value of those recorded failure/hazard limits, for example η j = max η X , η Y , η z , . . . .In the latter, target dynamic system survival proba- bility P( ) being defined as a function of parameter ; notice that P ≡ P(1) being derived from Eq. (1).Now it is possible to evaluate the target system's probability of non-exceedance (survival) P( ) directly Because dependence between neighboring R j may not always be negligible, condition- ing level k = 1 being imposed after one-step statistical memory approximation for 2 ≤ j ≤ N (conditioning level k = 2 ).Equation (4)'s approximation may also be fur- ther assessed where 3 ≤ j ≤ N (conditioning level k = 3 ) etc. Statistical correlations between close local maxima are better captured by the latter class of approximations p k ( ) on the conditioning level k and independent of j .Now, the modified (4-parameter) Weibull approach may be used to extrapolate non-exceedance (survival) probability (Gaidai et al. 2023a(Gaidai et al. , 2010a(Gaidai et al. , 2010b(Gaidai et al. , 2018(Gaidai et al. , 2016(Gaidai et al. , 2022a(Gaidai et al. , 2022b;;Balakrishna et al. 2022) Well-known mean up-crossing rate function can be used to express exceedance risk/ probability (Gaidai et al. 2022c(Gaidai et al. , 2022d, 2022e, 2022f;, 2022e, 2022f;Gaidai and Xing 1994;Stanisic et al. 2018;Zhang et al. 2018Zhang et al. , 2019;;Xu et al. 2019), given by Eq. ( 5).When the conditioning parameter k being large enough, one has (3) 6) for k = 1 representing mean up-crossing rate function For the previously presented non-dimensional synthetic vector R(t) , being put together, using scaled MDOF dynamic system's combined vector X η X , Y η Y , Z η Z , . . . .Next, ν + ( ) being the mean up-crossing rate function of the non-dimensional system failure/ hazard level .Although system stationarity has been presumed, the recommended technique could work well in non-stationary situations.Given an in-situ scatter diagram with m = 1, .., M environmental in-situ sea states, each of which has a probability of q m with M m=1 q m = 1 , with a matching long-term equation the same function as in Eq. ( 7), with p k ( , m) standing for a particular short-term ambi- ent sea condition, having index m.For high values of , the probability-related p k ( ) functions, shown above, approaching non-dimensional target level 1 as they are often quite regular in the PDF tail ( ≥ 0 ).PDF tail frequently exhibits behavior, resembling exp −(a + b) c + d , where 4 parameters a, b, c, d being constants fitted to the proper PDF tail cut-on 0 value.The sequential quadratic programming (SQP) approach, provided by Numerical Library from the Numerical Algorithms Group (NAG) (Gaidai et al. 2016) may be used to find optimal parameter values for the 4 parameters/variables a, b, c, d (Gaidai et al. 2022g, 2022h, 2022i, 2023c(Gaidai et al. 2022g, 2022h, 2022i, , 2022j, 2023d;;Gaidai andXing 2022a, 2022b).

Results
In this section we present validation of Gaidai reliability method, using both synthetic as well as onboard measured datasets.

Synthetic example
In this sub-section synthetic example being presented, as the exact analytical solutions for underlying probability distributions being known in advance.That will allow the authors to compare statistical methods and cross-validate them versus the exact analytical predictions.Let one consider the 3.65-dayly maximum windspeeds process X(t) simulated during given time lapse [0, T ] .The underlying normalised non-dimensional stochastic processes U (t) has been modelled as a Gaussian stationary process with 0 mean value and a standard deviation equal to 1 Hence, it was assumed that the U (t) mean zero up-crossing rates satisfies the equality ν + U (0) = 10 3 /T , with T = 1 year, the latter assumption is common in offshore wind engineering (Gaidai et al. 2022e).For this numerical example, the data record had been chosen to contain 10 4 data points, which is (7) equivalent to 100 years record, since windspeed maxima process X(t) has 365/3.65=10 2 data points per year.
The underlying windspeed process U (t) results in 3.65 days maximum analytical CDF (Cumulative Density distribution) distribution F 3d X (x) = exp −qexp − x 2 2 for the windspeed three days maxima process X 3d (t) , with q = 10 .There are 3 Archimedian copulas being in common use: Clayton, Frank and Gumbel-Haugaard.First, the Gumbel-Haugaard copula dependence copula G(u, v) between the two marginal peak wind- speed variables X 3d (t) and identically distributed correlated process Y 3d (t) being considered with parameter m = 1/ √ 1 − R corr being related to the correlation coefficient R corr between two processes X 3d (t) and Y 3d (t); in this section R corr was set to 0.5 for sim- plicity.As the underlying windspeed processes X(t) , Y (t) are identically distributed stationary Gaussian processes, and Gumbel-Haugaard copula has been relatively easy to fit for the 4-parameter bivariate (2D) Weibull method.It is expected that it will be no significant difference between Gaidai reliability and 4-parameter bivariate Weibull method predictions, in terms of predicted decimal logarithm probability level for a given response level (windspeed in this case) of interest.For this section, the following response levels have been chosen: x = 6 , y = 5.2 .2D extreme value distribution of local peak's event data was Figure 3 presents Monte Carlo simulated timeseries, along with corresponding 4-parameter bivariate Weibull method bivariate contour and Gaidai reliability method's (10) 3 Simulated timeseries generated using Gumbel-Haugaard copula prediction based on − → R vector, red star indicated target probability level, see Fig. 4. As expected, the agreement between both methods was quite good.Second, the Clayton copula C(u, v) being an asymmetric Archimedean copula, has been applied in an analo- gous way, instead of the Gumbel-Haugaard copula Clayton copula is less convenient to fit by 4-parameter bivariate Weibull method as it is not in its copula library.Therefore, in this case, 4-parameter bivariate Weibull method is expected to perform less accurately than the Gaidai reliability method, as obviously, Gaidai reliability does not have any copula approximation assumptions.All numerical parameters have been unchanged, with respect to the previous Gumbel-Haugaard copula case.Series of independent synthetic tests found that the Gaidai reliability method predicts the exceedance probability's decimal logarithm level more accurately than 4-parameter bivariate Weibull method.For the numerical setup described above, it was found that, on average, Gaidai reliability performs about 15% more accurately than 4-parameter 2D Weibull method.However, synthetic data was based on the underlying Gaussian process and Archimedian copulas.Gaidai reliability advantages may be much more pronounced in the case of real measured non-Gaussian, cross-correlated by non-Archimedian copula data.Last but not least, the computational effort of 4-parameter bivariate Weibull method being more significant than the Gaidai reliability method for any given bivariate failure limit, as 4-parameter bivariate Weibull method performs twodimensional surface interpolation; the latter is significant when analysing large data sets, as in the next section.

Cargo ship onboard analysis
This sub-section illustrates the Gaidai reliability method's efficiency by application to container vessel onboard measured ship panel stresses data set.It is known that container ship panel dynamic stresses due to container ship hull vibrations constitute a highly nonlinear multi-dimensional and cross-correlated dynamic system that is challenging to analyse.Moreover, this vessel dynamic system reliability study is paramount for container vessels crossing Atlantics and experiencing extreme weather conditions.( 12) Even though these two ships may not have been designed and approved according to best practise resulting in substandard collapse strength compared to other similar container ships, both ships broke due to hull girder overloading.MSC Napoli broke in the engine room bulkhead, and MOL Comfort broke amidships.Such severe accidents as mentioned in the above 2 cases, have been followed up by thorough investigations.Container 4400TEU Panamax ship has been equipped with motion sensors during its whole trans-Atlantic voyages instrumented by DNV (DNV-RP-H103 2011).Further explanation of the setup of the measuring systems is given in Storhaug and Moe (2007).Container vessel has been instrumented using standard hull monitoring systems following Figure 5 presents 4400TEU mid-onboard strain sensors placement along with observed crack positions.Sensor placement was done according to DNV container vessel rules and regulations.Regarding three mid-ship censor locations, indicated by circles in Fig. 7, only upper deck (upper circle) measurements were utilised in this study.Similarly, sensors were placed aft of the vessel, resulting in measured stresses in the longitudinal direction on a flat bar below the upper deck.In this section following abbreviations will be used: P-port, S-starboard side DL-deck longitudinal In order to unify both measured timeseries X, Y the following scaling has been per- formed according to Eq. ( 11), making both two responses non-dimensional and having the same failure/hazard limit equal to 1. Next, all local maxima from both two measured timeseries were merged into one single timeseries by keeping them in time non-decreasing order: − → R = (max{X 1 , Y 1 }, . . ., max{X N , Y N }) .In order to unify both measured timeseries X, Y the following scaling was performed Note that Fig. 6 refers only to a chosen stress couple (bivariate failure/hazard point).Thus, when changing the couple, the figure should be re-drawn.
Figure 7 left presents 4-parameter bivariate Weibull method bivariate contours for DLP and DLS ship panel stresses, (Gaidai et al. 2023g, 2023h, 2023i, 2023j, 2023k, 2023l;Gaidai and Xing 2023;Xing and Gaidai 2023).It is seen from Fig. 7 left that 4-parameter bivariate Weibull method fits different Gumbel copula to the measured data, and there is an inherent error due to a particular copula choice.For more details on the 4-parameter bivariate Weibull method, see (Zhang et al. 2023;Liu et al. 2023;Gaidai et al. 2023m, 2023n, 2023o;Yakimov et al. 2023a).For more details on the Gaidai reliability method, (Sun et al. 2023;Gaidai et al. 2023p;Yakimov et al. 2023b).Bivariate failure point η X = 120 MPa, η Y = 70 MPa was chosen as seen in Fig. 7 left, then synthetic vector − → R was obtained from DLP and DLS stresses, namely DLS stress record was scaled by multiplying with η X /η Y , see Fig. 6.The probability level log 10 p = −5 corresponding to this contour line was then compared with Gaidai reliability method estimate, see Fig. 7 right.It was found that 4-parameter bivariate Weibull method probability level estimate lied well within 95% CI (Confidence Interval), predicted by the Gaidai reliability method.As indicated with circles in Fig. 6, other stress couples were studied, and again, the Gaidai reliability method was well verified by 4-parameter bivariate Weibull method.
Gaidai reliability handles the system as a black box with an infinite number of random parameters, as this method does not analyse any random parameters.The number of dimensions is also infinite, as all of them are merged into 1D vector R(t) .The ship panel stress response spectrum at almost all sea states contains two peaks (bimodal spectrum).WF (wave-frequency) part is related to the wave-induced responses, and HF (high-frequency) part may be caused by the transient load or resonant vibrations, known as whipping/springing.Future studies aiming at the whipping phenomenon should analyse the entire stress response process as in this work and separately address low and high-pass filtered responses.
The novel methodology is of general purpose and suitable for any jointly stationary dynamic system.If such measurements are available, various failure mechanisms, such as fatigue failure, corrosion-type fracture failure or other ultimate collapse mechanisms related to the coupled deck system in the vessel, may be added as new dimensions to the dynamic system of interest.In the long run, the study aims to develop a robust code for ship routing design concerning fatigue accumulation and extreme loading.This is (13) an essential task, and the proposed reliability approach may serve as the basis for future work.It is recommended that more studies are needed on the accurate stochastic models for the wave environments in different ocean regions.The possible results will significantly benefit modern cargo ship design.

Discussion
This study has presented novel reliability methodology, validated it using cargo ship onboard dataset.For complex dynamic system options like direct Monte Carlo simulation approach, or extensive laboratory experiments, or in the best-case onboard measurements are mostly out of practical reach.Existing reliability methods, used in naval architecture are mostly limited to 2D systems, while suggested Gaidai reliability methodology does not have any limitation on the number of system dimensions.Classic reliability methods, dealing with timeseries, do not have the advantage of dealing efficiently with systems possessing high dimensionality and cross-correlation between different system responses.The essential advantage of the introduced methodology is its ability to study the reliability of high-dimensional dynamic systems.The method presented in this paper has been previously validated by application to a wide range of simulation models, but for only one-dimensional system responses and, in general, very accurate predictions were obtained.This study aimed to develop a general-purpose further, robust and simple-to-use multi-dimensional reliability method.When speaking of onboard measurements, those datasets are rarely available for public research presentation, due to confidentiality issues.This study, however has managed to utilize real, extensive, raw measured onboard dataset.This research intends to contribute to the development of new DNV safety rules and standards for contemporary maritime industry, thus potentially influencing future cargo transport insurance policy, as well as cargo vessel operational policies.

Conclusions
This study analysed both synthetic windspeed data set and ship dynamic response timeseries measured onboard during over 70 transatlantic voyages of 4400TEU Panamax container vessel in 2010's years.The theoretical reasoning behind the proposed method is given in detail.Note that using direct measurement or Monte Carlo simulation for dynamic system reliability analysis is attractive.However, dynamic system complexity and its high dimensionality require developing novel robust and accurate techniques that can deal with a limited data set at hand, utilising available data as efficiently as possible.The novel method was validated versus the 4-parameter bivariate Weibull method bivariate method using synthetic, analytical data and actual onboard measured ship panel stress data.Finally, the suggested methodology can be used in various engineering areas of applications.The presented naval architecture example does not limit areas of new method applicability.
For future work, the authors will include other container vessel types of a larger size and have more sensor measurements at other critical deck panel locations.Particular focus will be laid on larger TEU vessels whipping responses, simultaneously recorded at different vessel hull locations.When speaking about advocated reliability method

Table 1
Main particulars of 4400TEU vessel