Gaidai multivariate risk assessment method for cargo ship dynamics

ABSTRACT

Contemporary cargo vessel transport constitutes a vital part of the global trade and economy. Thus, it is of primary engineering and design importance to further develop more efficient novel risk analysis and risk assessment methods for large cargo ships. Classical risk assessment and risk analysis methodsmay not always have advantages of efficiently dealing with dynamic system’s high dimensionality along with cross-correlations between various critical dimensions. Current study benchmarks state-of-art Gaidai spatiotemporal structural risk assessment method, suitable for multivariate structural dynamic systems versus a well-established bivariate risk assessment method. As an example, for this risk assessment study, an operating container vessel has been selected. The vessel hull had been subjected to large deck panel stresses and accelerations during cross-Atlantic sailing. Risks of losing cargo containers are caused by extreme dynamics being one of the primary concerns for vessel cargo transport. Gaidai multivariate risk assessment methodology, advocated here, opens up new possibilities of predicting efficiently, yet simply potential structural damage/failure risks for nonstationary, nonlinear, multivariate cargo ship dynamic systems, as a whole. Note that advocated novel multivariate risk assessment method may be applied to a wide range of complex marine and offshore engineering systems.

KEYWORDS:

• Container ship
• trade
• AI
• ship motions
• dynamic system
• transportation

1. Introduction

Extreme value predictions problems in contemporary engineering being challenging and frequent, especially in case when measured data are scarce (Abioye et al., 2021; Chen et al., 2021; Drummen et al., 2008; Dulebenets, 2022; Elmi et al., 2022; European Maritime Safety Agency EMSA, 2021; ISSC, 2009; Jovanović et al., 2022; Parunov et al., ; Y. Wang & Wang, 2021). Thus, from practical design standpoint, it is essential to develop novel, effective, accurate risk assessment approaches. As container vessels play undeniably significant role within commercial marine transport industry, this study primarily has been concerned with cargo vessel’s safety and risk assessment. MSC Napoli, post-Panamax container vessel, had experienced engine room bulkhead’s damage in Jan 2007. Evaluations of MSC Napoli’s hull girder risk and risk assessment at the time of the accident have been presented by (Parunov et al.,). Hull girder ultimate strength was representative for vessel’s structural capacity; loads included still water vessel hull bending moments, vertical wave-induced bending moments, whipping-induced bending moments. Load combinations of wave bending moments and whipping have been taken into account, first-order risk assessment approach had been utilized to assess failure risks (Parunov et al.,). Another post-Panamax container vessel, MOL Comfort, had been damaged in Jun 2013 (Bureau Veritas, 2012; Ellermann, 2008; Falzarano et al., 2012; Gaidai, Yakimov, et al., 2023; Rahim, 2019; Su, 2012). Both cargo vessels have been damaged due to overloading of the main deck panel’s girder, resulting in a decreased collapse strength, compared to other comparable vessels. These significant failure events, involving cargo ships, have seriously rattled maritime industry, calling for thorough investigations. Next, one may mention human-related error risk assessment aspects, typical within maritime cargo transportation industry (Bicen et al., 2021). A wide range of intricate factors, such as insufficient experience, training, system monitoring, poor communication, or a failure to draw lessons from previous mistakes, can lead to human error.

This study’s major research object has been 4400 TEU (Twenty-foot Equivalent Unit) Panamax container vessel, 245-meter long (ISSC, 2009). Motion sensors have been installed on this 4400 TEU Panamax container vessel, during trans-Atlantic voyages during 2010–2020 years (Abioye et al., 2021; Drummen et al., 2008; Dulebenets, 2022; Elmi et al., 2022). Captain had made additional measures, in order to avoid worst incoming storms. 4400 TEU dynamic vessel panel stress sensor’s placements had been carried out following latest DNV (Det Norske Veritas) container vessel’s safety regulations and standards (Chen et al., 2021; Ellermann, 2008; European Maritime Safety Agency EMSA, 2021; Falzarano et al., 2012; Jovanović et al., 2022; Parunov et al., ; Rahim, 2019; Su, 2012). On-board data’s sampling time-step has been $\mathit{dt}=0,025$ s, monitored continuously through complete voyage’s duration. For more details on vessel’s parameters, see e.g. (Bicen et al., 2021; Bureau Veritas, 2012; DNV-RP-H103, 2011; Gaidai, Yakimov, et al., 2023). On-board measured dataset had not been post-processed/filtered, hence raw dataset had been analyzed; all available voyage-datasets have been coalesced into a single 1D (1-Dimensional), jointly representative timeseries. Evaluating extreme vessel accelerations and deck panel’s stresses along with associated low risks/probabilities being challenging task, only limited progress had been made within the past decade (Andersen & Jensen, 2014; Gaidai et al., 2010a; Gaidai, Yakimov, et al., 2023; Storhaug & Moe, 2007; Sun, Gaidai, Xing, et al., 2023; Tomic et al., 2018). Numerical MC (Monte Carlo) simulations, being based on underlying model’s assumptions, often being unable to capture all necessary wave-induced accelerations and vessel deck panel’s stresses.

Figure 1. Example of cargo-loaded TEU cargo ship.

When assessing design loads given current vessel regulations, vessel hulls being mostly modelled as rigid. Cargo vessels being in fact flexible, with waves causing their structural hull vibrations. Environmental loads, along with accumulated fatigue damage, are being exacerbated due to vessel hull vibrations. Full-scale measurements for Capesize bulk type carriers have provided evidences of key risk contributions, but impacts on container overall’s overall risk assessment had not been yet comprehensively studied. 294-m-long cargo vessel sailing in North Atlantic had been used in this study in order to assess areal deck panel stresses, along with vessel’s accelerations. Wave radar along with wind data had been utilized to complement stress’s records amidships, those measurements have been evaluated aboard, in real time, thanks to hull onboard monitoring system, giving captain updated live information. Axial forces, warping, horizontal as well as vertical bending moments all are contributing to increased stress levels, along with fatigue damage. Differences between overall fatigue damage, compared to wave frequency-related damage being typically utilized to compute fatigue damage, caused by wave-induced hull vibrations. Wave-induced hull vibrational effects are depending on vessel’s size, its eigen-periods, built materials, initial cracks. To decrease maintenance’s downtime throughout whole course of ship’s service lifetime, wave-induced vibrational effects have to be taken into account at the design verification stage for any large cargo ship. On-board measured dataset being vital to correctly assess structural vessel’s risk assessment, as well as calibrate software numerical tools, contribute to vessel design rules/regulations.

Constructed in the year 2003, 4400 TEU vessel has followed a contemporary design. DNV has assessed on-board stresses and acceleration records from a 4400 TEU container vessel having comparable size, age as a part of AOG (Active Operator Guidance) DNV research project. Full-scale type measurements, recorded in North Pacific exhibited contribution up to 50% depending on season, in contrast to about 35% according to head sea model tests (Drummen et al., 2008), Figure 2.

Figure 2. Sensor, being glued to the deck panel’s stiffener web (Storhaug & Moe, 2007).

4400 TEU vessel full scale on-board measurements, along with wave-induced vibrational analysis tools have been provided by the DNV.

Table 1. 4400 TEU vessel main particulars.

Length overall, LOA [m]	294.0
Length between perpendiculars, LPP [m]	281.0
Breadth, B [m]	32.3
Depth, D [m]	21.5
Draft’s design, T [m]	10.8
Draft’s scantlings, T [m]	10.8
Deadweight at design, dwt [Megatonne]	0.5
Service speed at design draft [knot]	23.0
Block coefficient, CB [-]	0.7
Power [MW]	37.3

Figure 1 presents 4400 TEU cargo vessel, exploited in a trade between Germany, UK, France and Canada. Table 1 presents brief summary of vessel’s essentials. Vessel received notations: +1A1, EO, Container Carrier, NAUTICUS (New-building), W1-OC, ICE-C, DG-P, constructed according to DNV class. Vessel had typical cross-sectional shape, however running with somewhat lower draft, compared to other ships of one’s size. HT32 steel was utilized, with HT36 above top deck, hatch coaming plating was 65 mm in thickness, shear strake along with top deck were built of 60 mm thick steel plates. Sensors SW, IS, indicated in Figure 5 have been intended to be placed within stiffener web’s neutral axis, at stiffener’s momentary zero point, to monitor bending caused by lateral pressures.

In the next Section, state-of-art Gaidai multivariate risk assessment methodology to be introduced in some detail, following by practical application of Gaidai multivariate risk assessment method to 4400 TEU cargo ship’s on-board measured data.

2. Gaidai multivariate risk assessment methodology

Assessing structural vessel system’s safety, utilizing conventional design/engineering risk assessment methods being often quite challenging (Abioye et al., 2021; Chen et al., 2021; Drummen et al., 2008; Dulebenets, 2022; Elmi et al., 2022; European Maritime Safety Agency EMSA, 2021; Falzarano et al., 2012; Jovanović et al., 2022; Parunov et al., ; Rahim, 2019; Y. Wang & Wang, 2021). Risk assessment of complex multivariate structural vessel systems may be directly assessed either by conducting extensive direct numerical MC simulations or by possessing sufficiently representative simulated/measured dataset. Most of the dynamic nonstationary engineering marine systems being often too complex for be tackled by both direct MC computational, as well as experimental lab/onboard methods. Authors have recently developed state-of a-t multivariate risk assessment approach, suitable for structural vessel hull systems in an effort to further reduce computational/measurement costs.

Ocean waves have been typically thought as following quasi-stationary, homogenous, piecewise-ergodic stochastic process, having 3-h storms stationary ‘windows’. Think of an offshore engineering or naval structure having large number of degrees of freedom, subjected to random external environmental in situ loadings, that are temporally quasi-stationary, e.g. environmental in situ waves and wind. Another option being to regard dynamic process as being reliant on a number of covariates, whose dynamic change over time may be regarded a quasi-ergodic process. Dynamic vessel system’s vector being now composed/assembled of either measured or numerically MC simulated$\left(\mathit{X}\left(\mathit{t}\right),\mathit{Y}\left(\mathit{t}\right),\mathit{Z}\left(\mathit{t}\right),\dots \right)$ vessel hull system’s critical/key components $\mathit{X}\left(\mathit{t}\right),\mathit{Y}\left(\mathit{t}\right),\mathit{Z}\left(\mathit{t}\right),\dots$ synchronously recorded though sufficient (representative) timelapse $\left(0,\mathit{T}\right)$. Dynamic vessel’s system 1D key/critical component’s global maxima, being measured through whole (representative) timelapse $\left(0,\mathit{T}\right)$ and denoted here as $X_T^{\mathrm{m}\mathrm{a}\mathrm{x}}={max}_{0\le \mathit{t}\le \mathit{T}}\mathit{X}\left(\mathit{t}\right)$, $Y_T^{\mathrm{m}\mathrm{a}\mathrm{x}}={max}_{0\le \mathit{t}\le \mathit{T}}\mathit{Y}\left(\mathit{t}\right)$, $Z_T^{\mathrm{m}\mathrm{a}\mathrm{x}}={max}_{0\le \mathit{t}\le \mathit{T}}\mathit{Z}\left(\mathit{t}\right),\dots$, respectively. With large-enough $\mathit{T}$ value one means large value, compared to dynamic vessel hull system’s relaxation, as well as auto-correlation temporal scales (Andersen & Jensen, 2014; Balakrishna et al., 2022; Gaidai et al., 2010a, 2010b, 2016, 2018, 2022b; Gaidai, Xu, et al., 2022; Gaidai, Yakimov, et al., 2023, 2023; Gaidai, Yan, et al., 2022; Storhaug & Moe, 2007; Sun, Gaidai, Xing, et al., 2023). Let $X_1,\dots ,X_{N_X}$ be temporally subsequent dynamic vessel hull system critical/key component’s process’s $\mathit{X}\left(\mathit{t}\right)$local maxima, measured/recorded at discrete time-instants $t_1^X<\dots <t_{N_X}^{\mathit{X}}$ within the observational time period $\left(0,\mathit{T}\right)$. System’s critical/key components $\mathit{X}\left(\mathit{t}\right),\mathit{Y}\left(\mathit{t}\right),\mathit{Z}\left(\mathit{t}\right),\dots$ being either vessel system’s key/critical loads, or vessel hull’s dynamic responses, on their own. Local maxima for the rest of MDOF dynamic vessel system’s components, e.g. $\mathit{Y}\left(\mathit{t}\right),\mathit{Z}\left(\mathit{t}\right),\dots$, being designated as $Y_1,\dots ,Y_{N_Y};$ $Z_1,\dots ,Z_{N_Z}$ etc. Ship system’s key/critical component’s local maxima had been assumed being positive. Assessing target damage/hazard/failure risk/probability for a dynamic vessel hull system being now primary objective.

(1) $1-\mathit{P}=\mathrm{Prob}(X_T^{\mathrm{m}\mathrm{a}\mathrm{x}}>{\eta }_X\cup Y_T^{\mathrm{m}\mathrm{a}\mathrm{x}}>{\eta }_Y\cup Z_T^{\mathrm{m}\mathrm{a}\mathrm{x}}>{\eta }_Z\cup \dots )$(1)

with

(2) $\mathit{P}=\underset{\left(0,0,0,\dots \right)}{\overset{\left({\eta }_X,{\eta }_Y,{\eta }_Z,\dots \right)}{\int \int \int }}{p}_{X_T^{\mathrm{m}\mathrm{a}\mathrm{x}},Y_T^{\mathrm{m}\mathrm{a}\mathrm{x}},Z_T^{\mathrm{m}\mathrm{a}\mathrm{x}},\dots }\left(X_T^{\mathrm{m}\mathrm{a}\mathrm{x}},Y_T^{\mathrm{m}\mathrm{a}\mathrm{x}},Z_T^{\mathrm{m}\mathrm{a}\mathrm{x}},\dots \right)\mathit{d}{X}_T^{\mathrm{m}\mathrm{a}\mathrm{x}}\mathit{d}{Y}_T^{\mathrm{m}\mathrm{a}\mathrm{x}}\mathit{d}{Z}_T^{\mathrm{m}\mathrm{a}\mathrm{x}}\dots$(2)

being target vessel system’s survival (non-exceedance) chances/probability for critical/key values of dynamic vessel hull system’s critical components ${\eta }_X$, ${\eta }_Y$, ${\eta }_Z$,… with $\cup$ denoting logical unity operator «or»; $p_{X_T^{\mathrm{m}\mathrm{a}\mathrm{x}},Y_T^{\mathrm{m}\mathrm{a}\mathrm{x}},Z_T^{\mathrm{m}\mathrm{a}\mathrm{x}},\dots }$ representing the joint PDF (Probability Density Function) of global vessel system’s key/critical component’s maxima, within timelapse $\left(0,\mathit{T}\right)$. Because of high dimensionality possessed by most dynamic engineering systems, along with underlying dataset limitations, in practice being is not possible to straightforwardly assess later vessel hull joint system’s PDF, $p_{X_T^{\mathrm{m}\mathrm{a}\mathrm{x}},Y_T^{\mathrm{m}\mathrm{a}\mathrm{x}},Z_T^{\mathrm{m}\mathrm{a}\mathrm{x}},\dots }$. Dynamic vessel system being regarded to have been instantaneously damaged/failed (or have entered state of hazard) at the temporal instants, when either $\mathit{X}\left(\mathit{t}\right)$ exceeding ${\eta }_X$, or $\mathit{Y}\left(\mathit{t}\right)$ exceeding ${\eta }_Y$, or $\mathit{Z}\left(\mathit{t}\right)$ exceeding ${\eta }_Z$, etc., fixed damage/hazard/failure levels ${\eta }_X$, ${\eta }_Y$, ${\eta }_Z$,… being potentially different for each single dynamic vessel system’s 1D component. $X_{N_X}^{\max }=max\{X_j;\mathit{j}=1,\dots ,N_X\}=X_T^{\mathrm{m}\mathrm{a}\mathrm{x}}$, $Y_{N_Y}^{\max }=\max \{Y_j;\mathit{j}=1,\dots ,N_Y\}=Y_T^{\mathrm{m}\mathrm{a}\mathrm{x}}$, $Z_{N_z}^{\max }=max\{Z_j;\mathit{j}=1,\dots ,N_Z\}=Z_T^{\mathrm{m}\mathrm{a}\mathrm{x}}$, and so forth. Temporal instants corresponding to local maxima $\left[t_1^X<\dots <t_{N_X}^{\mathit{X}};\mathit{\hspace{0.17em}}{t}_1^Y<\dots <t_{N_Y}^{\mathit{Y}};\mathit{\hspace{0.17em}}{t}_1^Z<\dots <t_{N_Z}^{\mathit{Z}}\right]$ being now coalesced into one single temporal synthetic system’s vector $t_1\le \dots \le t_N$ preserving temporally increasing order, with $t_N=\max \{t_{N_X}^{\mathit{X}},t_{N_Y}^{\mathit{Y}},t_{N_Z}^{\mathit{Z}},\dots \}$, $\mathit{N}\le N_X+N_Y+N_Z+\dots$ Occurrence times of vessel system’s key/critical component’s local maxima have been denoted here by $t_j$. MDOF risk/damage/limit/hazard vector $\left({\eta }_X,\mathit{\hspace{0.17em}}{\eta }_Y,\mathit{\hspace{0.17em}}{\eta }_Z,\dots \right)$, have been introduced, consisting of hazard/damage/failure limits of dynamic vessel hull system’s key/critical components$\mathit{X},\mathit{Y},\mathit{Z},\dots$. 1D dynamic system’s key load/response critical/key 1D component’s local maxima being now combined into single synthetic vector $\mathit{R}\left(\mathit{t}\right)\equiv \vec{R}=\left(R_1,R_2,\dots ,R_N\right)$, following merged system’s temporal vector $t_1\le \dots \le t_N$. In order to reduce limit/damage/hazard values for all load/response critical component, scaling parameter $0<\mathit{\lambda }\le 1$ to be introduced. MDOF limit/damage/hazard vector $\left({\eta }_X^{\lambda },{\eta }_Y^{\lambda },\mathit{\hspace{0.17em}}{\eta }_z^{\lambda },\dots \right)$ to be now introduced, having ${\eta }_X^{\lambda }\equiv \mathit{\lambda }\cdot {\eta }_X$, $\equiv \mathit{\lambda }\cdot {\eta }_Y$, ${\eta }_z^{\lambda }\equiv \mathit{\lambda }\cdot {\eta }_Z$, … The latter introduces the unified vessel hull system’s limit/hazard/damage/risk vector $\left({\eta }_1^{\lambda },\dots ,{\eta }_N^{\lambda }\right)$, having each critical component ${\eta }_j^{\lambda }$ being equal to either ${\eta }_X^{\lambda }$, or ${\eta }_Y^{\lambda }$or ${\eta }_z^{\lambda }$, … , depending on which specific vessel system’s critical/key component has supplied given local maxima, recorded at given temporal moment $t_j$; in case of simultaneous occurrence of few vessel hull system critical/key component’s local maxima at specific temporal instant $t_j$, values of ${\eta }_j^{\lambda }$ have been chosen here to be equal maximum value of those recorded failure/damage/hazard limits, e.g. ${\eta }_j^{\lambda }=\max \left\{{\eta }_X^{\lambda },\mathit{\hspace{0.17em}}{\eta }_Y^{\lambda },\mathit{\hspace{0.17em}}{\eta }_z^{\lambda },\dots \right\}$. In the above, target dynamic vessel hull system’s survival chances/probability $\mathit{P}\left(\mathit{\lambda }\right)$ being defined as a regular function of the scaling parameter $\mathit{\lambda }$; with $\mathit{P}\equiv \mathit{P}\left(1\right)$ following Eq. (1). Now it is well possible to assess target vessel hull system’s non-exceedance (survival) probability $\mathit{P}\left(\mathit{\lambda }\right)$ directly

(3) $\mathit{P}\left(\mathit{\lambda }\right)\mathit{\hspace{0.17em}}=\mathrm{Prob}\left\{R_N\le {\eta }_N^{\lambda },\dots ,R_1\le {\eta }_1^{\lambda }\right\}=\mathrm{Prob}\{R_N\le {\eta }_N^{\lambda }|R_{\mathit{N}-1}\le {\eta }_{N-1}^{\lambda },\dots ,R_1\le {\eta }_1^{\lambda }\}\cdot \mathrm{Prob}\left\{R_{\mathit{N}-1}\le {\eta }_{N-1}^{\lambda },\dots ,R_1\le {\eta }_1^{\lambda }\right\}=\prod _{\mathit{j}=2}^N\mathrm{Prob}\{R_j\le {\eta }_j^{\lambda }|R_{\mathit{j}-1}\le {\eta }_{1j-}^{\lambda },\dots ,R_1\le {\eta }_1^{\lambda }\}\cdot \mathrm{Prob}\left(R_1\le {\eta }_1^{\lambda }\right)$(3)

Due to inter-dependencies between the neighboring $R_j$ are not always be negligible, level of conditioning $\mathit{k}=1$ to be imposed after 1-step statistical/risk assessment memory’s approximation

(4) $\mathrm{Prob}\{R_j\le {\eta }_j^{\lambda }|R_{\mathit{j}-1}\le {\eta }_{j-1}^{\lambda },\dots ,R_1\le {\eta }_1^{\lambda }\}\approx \mathrm{Prob}\{R_j\le {\eta }_j^{\lambda }|R_{\mathit{j}-1}\le {\eta }_{j-1}^{\lambda }\}$(4)

for $2\le \mathit{j}\le \mathit{N}$ (level of conditioning $\mathit{k}=2$). Eq. (4)‘s approximation to be further assessed

(5) $\mathrm{Prob}\{R_j\le {\eta }_j^{\lambda }|R_{\mathit{j}-1}\le {\eta }_{j-1}^{\lambda },\dots ,R_1\le {\eta }_1^{\lambda }\}\approx \mathrm{Prob}\{R_j\le {\eta }_j^{\lambda }|R_{\mathit{j}-1}\le {\eta }_{j-1}^{\lambda },R_{\mathit{j}-2}\le {\eta }_{j-2}^{\lambda }\}$(5)

where $3\le \mathit{j}\le \mathit{N}$ (level of conditioning $\mathit{k}=3$) etc. Statistical inter-correlations due to clustered/grouped local maxima may be better captured by the latter approximations’ series$p_k\left(\mathit{\lambda }\right):=\mathrm{Prob}\{R_j>{\eta }_j^{\lambda }|R_{\mathit{j}-1}\le {\eta }_{j-1}^{\lambda },R_{\mathit{j}-\mathit{k}+1}\le {\eta }_{j-k+1}^{\lambda }\}$with $\mathit{j}\ge \mathit{k}$ being dependent only on level of conditioning $\mathit{k}$ and being independent of $\mathit{j}$. Now, modified (4-parameter) Weibull’s approach to be utilized for survival (non-exceedance) probability tail extrapolation (Balakrishna et al., 2022; Gaidai et al., 2010a, 2010b, 2016, 2018, 2022b; Gaidai, Yakimov, et al., 2023; Gaidai, Yan, et al., 2022)

(6) $P_k\left(\mathit{\lambda }\right)\approx \exp (-\mathit{N}\cdot p_k\left(\mathit{\lambda }\right)),\mathit{k}\ge 1,\mathit{N}\gg \mathit{k}$(6)

MUR (Mean Up-crossing Rate) function may be utilized when expressing exceedance probability/risk (Gaidai, Fu, et al., 2022; Gaidai, Wang, et al., 2022; Gaidai & Xing, 2022c; Gaidai, Xu, et al., 2022; Gaidai et al., 2022a; Stanisic et al., 2018; X. Xu et al., 2019; J. Zhang et al., 2018, 2019), following Eq. (5). If conditioning parameter $\mathit{k}$ being large

(7) $\mathit{P}=\underset{\mathit{k}\to \mathrm{\infty }}{lim}{P}_k\left(1\right);\mathit{p}\left(\mathit{\lambda }\right)=\underset{\mathit{k}\to \mathrm{\infty }}{lim}{p}_k\left(\mathit{\lambda }\right)$(7)

Eq. (6) for $\mathit{k}=1$ representing MUR function

(8) $\mathit{P}\left(\mathit{\lambda }\right)\mathit{\hspace{0.17em}}\approx \exp (-{\nu }^{+}\left(\mathit{\lambda }\right)\mathit{T});{\nu }^{+}\left(\mathit{\lambda }\right)={\int }_0^{\mathrm{\infty }}\mathit{\zeta }{p}_{R\dot{R}}\left(\mathit{\lambda },\mathit{\zeta }\right)\mathit{d\zeta }$(8)

Non-dimensional synthetic vector that was previously introduced $\mathit{R}\left(\mathit{t}\right)$, being represented as scaled MDOF dynamic vessel system’s coalesced vector $\left(\frac{X}{{\eta }_X},\frac{Y}{{\eta }_Y},\frac{Z}{{\eta }_Z},\dots \right).$ Next, ${\nu }^{+}\left(\mathit{\lambda }\right)$ being MUR function of non-dimensional vessel system’s failure/damage/hazard level $\mathit{\lambda }$. Although vessel hull system’s quasi-stationarity had been pre-assumed, Gaidai multivariate risk assessment methodology being equally suitable for non-stationary systems. Provided in-situ scatter diagram, having $\mathit{m}=1,..,\mathit{M}$ environmental in-situ ocean/sea states, with each of which having probability/chance of occurrence of $q_m$ with ${\sum }_{m=1}^M{q}_m=1$, matching the long-term equality

(9) $p_k\left(\mathit{\lambda }\right)\equiv \sum _{\mathit{m}=1}^M{p}_k\left(\mathit{\lambda },\mathit{m}\right)q_m$(9)

with $p_k\left(\mathit{\lambda },\mathit{m}\right)$ matching particular short-term ambient ocean/sea condition, having index m. For high/extreme values of $\mathit{\lambda }$ (approaching non-dimensional target level $\mathit{\lambda }=1$) probability-related $p_k\left(\mathit{\lambda }\right)$ functions, introduced above, being often quite regular within PDF’s tail ($\mathit{\lambda }\ge {\lambda }_0$). PDF frequently exhibits tail behavior, often resembling $\exp \left\{-{\left(\mathit{a\lambda }+\mathit{b}\right)}^c+\mathit{d}\right\}$, with four parameters $\mathit{a},\mathit{b},\mathit{c},\mathit{d}$ being constants/parameters, fitted to roper PDF tail’s cut-on ${\lambda }_0$ value. SQP (Sequential Quadratic Programming) approach, created by NAG (Numerical Algorithms Group) library (Gaidai et al., 2016), may be well utilized in order to optimally fit parameters/variables $\mathit{a},\mathit{b},\mathit{c},\mathit{d}$, optimal parameter values, see (Gaidai & Xing, 2022a, 2022b; Gaidai, Wang, et al., 2022; Gaidai, Wu, et al., 2022; Gaidai, Xing, et al., 2023; Gaidai, Xu, et al., 2022, 2022; Gaidai, Yan, et al., 2023).

3. Results

This section presents Gaidai multivariate risk assessment method’s validation, utilizing synthetic, along with analogue measured on-board datasets.

3.1. Synthetic MC example

Current section presents synthetic MC example, when exact underlying PDF analytical solution being known. That allows statistical methods direct validation versus the exact analytical prediction. Let us consider 3.65-dayly maxima windspeeds process $\mathit{X}\left(\mathit{t}\right)$ that has been MC simulated, across the whole representative timelapse$\left[0,\mathit{T}\right]$. Underlying non-dimensionalized normalized stochastic/random processes$\mathit{U}\left(\mathit{t}\right)$ has been modelled here as Gaussian stationary process, having zero mean, unit standard deviation. Thus, it had been assumed that $\mathit{U}\left(\mathit{t}\right)$ MUR function satisfies equality${\nu }_U^{+}\left(0\right)={10}^3/\mathit{T}$, where $\mathit{T}=1$ year. The latter artificial assumption being not uncommon within wind and offshore engineering (Gaidai, Fu, et al., 2022). In this MC numerical example, the dataset contained ${10}^4$ datapoints, equivalent to a 100 years record, as windspeed’s maxima process $\mathit{X}\left(\mathit{t}\right)$ has 365/3.65=${10}^2$ datapoints per year. Underlying windspeed’s process $\mathit{U}\left(\mathit{t}\right)$ resulting in 3.65 days maxima CDF (Cumulative Distribution distribution) $F_X^{3d}\left(\mathit{x}\right)=\exp \left\{-\mathit{q}\exp \left(-\frac{x^2}{2}\right)\right\}$for windspeed’s 3-day maxima process $X^{3\mathit{d}}\left(\mathit{t}\right)$, with $\mathit{q}=10$. 3 Archimedean dependence copulas, being commonly in use: Clayton, Gumbel-Haugaard, Frank. 1^st, Gumbel-Haugaard’s copula $\mathit{G}\left(\mathit{u},\mathit{v}\right)$ between 2 peak windspeed’s marginal variables $X^{3\mathit{d}}\left(\mathit{t}\right)$ and identically distributed cross-correlated stochastic process $Y^{3\mathit{d}}\left(\mathit{t}\right)$ to be generated

(10) $\mathit{G}\left(\mathit{u},\mathit{v}\right)=\exp \left\{-{\left[{\left(-\log u\right)}^m+{\left(-\log v\right)}^m\right]}^{\frac{1}{m}}\right\}$(10)

having parameter $\mathit{m}=1/\sqrt{1-R_{\mathit{corr}}}$ expressed through correlation coefficient $R_{\mathit{corr}}$ between 2 processes $X^{3\mathit{d}}\left(\mathit{t}\right)$ and $Y^{3\mathit{d}}\left(\mathit{t}\right);$ in this section $R_{\mathit{corr}}$ has been set to 0.5. Underlying windspeed’s processes $\mathit{X}\left(\mathit{t}\right)$, $\mathit{Y}\left(\mathit{t}\right)$ being identically distributed and stationary Gaussian processes, Gumbel-Haugaard dependence copula had been quite easy to fit to 4-parameter 2D (bivariate) Weibull’s method. It was expected that there will be no any significant differences between the Gaidai multivariate risk assessment method and 4-parameter bivariate Weibull’s method forecasts, in terms of the predicted decimal log probability level, for a given critical/risk level of interest. In this section, following critical bivariate level had been chosen: $\mathit{x}=6$, $\mathit{y}=5.2$. 2D extreme value CDF corresponding to local peak events dataset being

(11) $H^{3\mathit{d}}\left(\mathit{x},\mathit{y}\right)=\exp \left\{-{\left[\mathit{q}\exp \left(-m\frac{x^2}{2}\right)+\mathit{q}\exp \left(-\mathit{m}\frac{y^2}{2}\right)\right]}^{\frac{1}{m}}\right\}$(11)

Figure 3 illustrates MC-simulated stochastic timeseries, along with 4-parameter 2D Weibull method contour, and Gaidai multivariate risk assessment method’s forecast, based on $\vec{R}$ vector, the star in Figure 4 indicates the target risk/probability level. As has been expected, agreement between two methods has been quite good. 2^nd, Clayton dependence copula $\mathit{C}\left(\mathit{u},\mathit{v}\right)$, being an asymmetric Archimedean copula, had been applied in a similar way, instead of the Gumbel-Haugaard’s dependence copula

Figure 3. MC Timeseries, generated utilizing Gumbel-Haugaard dependence copula.

Figure 4. Left: 4-parameter 2D Weibull’s contour. Right: matching Gaidai multivariate risk assessment method’s forecast, based on system’s synthetic $\vec{R}$ vector.

(12) $\mathit{C}\left(\mathit{u},\mathit{v}\right)=\max \left\{{\left[{\mathit{u}}^{-\mathit{m}}+v^{-\mathit{m}}\right]}^{-1\text{overm}},0\right\}$(12)

Clayton’s copula being less convenient for fitting by 4-parameter 2D Weibull’s method, as it was not within ‘copula library’. Thus, in the latter case, 4-parameter 2D Weibull method being expected to perform less accurate, than Gaidai multivariate risk assessment method, since Gaidai multivariate risk assessment does not rely on any dependence copula’s assumptions. All numerical parameters had been unchanged, wrt previous Gumbel-Haugaard’s copula case. Series of synthetic independent MC tests have shown that Gaidai multivariate risk assessment method predicted exceedance probability/risk decimal log-level even more accurately, than 4-parameter 2D Weibull’s method. For the numerical setup, described above, it has been found that Gaidai multivariate risk assessment method performed nearly 20% more accurately, than 4-parameter 2D Weibull’s method. Synthetic dataset has been based on underlying Gaussian process, as well as Archimedean copulas. Gaidai multivariate risk assessment methodology key advantages may be way more pronounced in case of real field measured non-Gaussian, inter-correlated by non-Archimedean dependence copula dataset. Computational effort, attributed to 4-parameter 2D Weibull’s method being way heavier, compared to Gaidai multivariate risk assessment methodology, for any given 2D failure/damage limit, as 4-parameter 2D Weibull’s method performing 2D surface’s interpolation. The latter being significant, if analysing large datasets, as in will be shown in the next section.

3.2. Cargo vessel on-board analysis

This section benchmarks Gaidai multivariate risk assessment method’s efficiency, proven by means of application to specific container ship on-board measured acceleration and deck panel stresses dataset. It has been known that container vessel panel’s dynamic stresses and accelerations, coupled to cargo vessel hull’s vibrations constituting highly nonlinear, multivariate, cross-correlated dynamic system, that is quite challenging to study. Cargo vessel dynamic system’s risk assessment study being important for large container ships, crossing Atlantic’s, experiencing potentially extreme ocean/weather conditions. MSC Napoli had been broken in vessel’s engine room’s bulkhead, while MOL Comfort had been broken amidships. Although these two ships may not have been designed/approved according to the best practise, recommended for collapse strength, compared to other similar container vessels, both vessels had been broken due to ship hull girder overloading. 4400 TEU Panamax container vessel had been instrumented with numerous motion on-board sensors during all trans-Atlantic voyages, DNV-instrumented (DNV-RP-H103, 2011). Further explanation of the measuring systems setup is given in (Storhaug & Moe, 2007). Container ship has been instrumented, using standard hull monitoring systems, in accordance with DNV hull monitoring standards/rules (DNV-RP-C205, 2010). Vessel has operated in trading activities between Germany, Belgium and Canada. Ship has not been built following DNV class but was transferred to DNV class, in the year 2004 receiving notation + 1A1, ICE-1A, EO, NAV-O, Container Carrier, TMON. The ship had been built in the year 1998, having pronounced bow flare, having somewhat old-fashioned cross-sectional vessel design. Two different ship panel stresses, located mid and aft of the ship have been chosen as critical components $\mathit{X},\mathit{Y},$ hence constituting practical example of 2D dynamic cargo ship dynamic system.

Figure 5 presents 4400 TEU strain sensors on-board placements, along with observed initial crack positions. Sensor’s placement has been done following DNV container ship’s regulations/rules. Only the upper vessel deck (see upper circle) sensor measurements have been utilised for this study. Analogously, sensors had been placed aft of the ship, resulting in on-board measured stresses following longitudinal direction on the flat bar, below upper deck. In this section, the following abbreviations to be utilized:

Figure 5. Locations of 4 amidship deck panel strain’s sensors, see Figure 2.

P – Port; S – Starboard

DL – Deck Longitudinal

DLP – Deck Longitudinal Port sensor (port side)

AFV – vessel’s forward acceleration

For unification of two measured timeseries $\mathit{X},\mathit{Y}$, scaling had been performed, following Eq. (11), making two key components non-dimensional with identical failure/damage/hazard limits, equal to 1. All local maxima, extracted from two measured timeseries have been merged/coalesced into one single timeseries, by keeping them in a temporally increasing order: $\vec{R}=\left(\max \left\{X_1,Y_1\right\},\dots ,\max \left\{X_N,Y_N\right\}\right)$. In order to unify two measured timeseries $\mathit{X},\mathit{Y}$, following scaling had been performed

(13) $\mathit{X}\to \frac{X}{{\eta }_X},\mathit{Y}\to \frac{Y}{{\eta }_Y}$(13)

with ${\eta }_X$, ${\eta }_Y$ being dimensional maxima, recorded over the whole observational timelapse, for dynamic vessel system’s key/critical components $\mathit{X},\mathit{Y}$, respectively.

Figure 6 refers only to a chosen stress-acceleration couple (2D failure/damage/hazard point).

Figure 6. Synthetic vector $\vec{R}$, derived from non-dimensionalized DLP ($\mathit{X}$) and AFV ($\mathit{Y}$) vessel dynamic components.

Figure 7 left presenting 4-parameter 2D Weibull’s method contours for DLP, AFV (Gaidai & Xing, 2023; Gaidai, Cao, Xing, et al., 2023; Gaidai, Cao, Xu, et al., 2023; Gaidai, Wang, et al., 2023; Gaidai, Xing, et al., 2023; Gaidai, Xu, et al., 2023, 2023; Xing & Gaidai, 2023). As seen from Figure 7 left, 4-parameter 2D Weibull’s method fits various Gumbel copula to the underlying measured dataset; hence, there is an inherent error, due to particular copula choice. For more details on 4-parameter 2D Weibull’s method, see (Gaidai, Xing, et al., 2023; Gaidai, Yan, et al., 2023; Gaidai et al., 2023a; Liu et al., 2023; Yakimov, Gaidai, Wang, & Wang, 2023; J. Zhang et al., 2023). For further details on Gaidai multivariate risk assessment method, [77]- (Gaidai, Yakimov, et al., 2023). 2D failure/hazard non-dimensional design point ${\eta }_1=0.79$, ${\eta }_2=0.55$ had been chosen, as seen in Figure 7 left; synthetic vector $\vec{R}$ has been derived from non-dimensionalized DLP and AFV, namely, AFV acceleration record has been scaled by means of multiplying with factor ${\eta }_1/{\eta }_2$, see Figure 6. Probability/risk threshold/level $\mathrm{lo}{\mathrm{g}}_{10}\mathit{p}=-4.5$ corresponding to this particular contour line had been then compared with Gaidai multivariate risk assessment method’s estimate, see Figure 7 right. It has been found that modified (4-parameter) 2D Weibull method risk/probability level forecast lied well within 95% CI, predicted by Gaidai multivariate risk assessment method, [87]

Figure 7. a) 4-parameter 2D Weibull’s method contours for non-dimensionalized DLP and AFV; b) corresponding Gaidai multivariate risk assessment forecast, dashed lines indicate 95% CI (confidence intervals).

The SODP (second-order difference plot) originated from the Poincare-type plot. SODP-type pllot reflects statistical patterns, residing within consecutive differences within the underlying timeseries dataset.

Figure 8 presents a second-order SODP-type plot, these kinds of plots may be used for dataset pattern recognition, QC (Quality Control), and comparison with other relevant similar datasets, for example using entropy AI (Artificial Intelligence) recognition approach (Gaidai et al., 2024; Gaidai et al., 2010a; Gaidai et al., 2010b; Gaidai et al., 2022; Gaidai et al., 2023b; Gaidai et al., 2023; Sun et al., 2023; Yakimov et al., 2023; Gaidai, Sheng, Cao, Zhu, et al., 2024).

Figure 8. Second (left) and 6^th order poincare difference plots.

Gaidai multivariate risk assessment method treats ship’s dynamic system as black box, having large number of covariates/random parameters NDOF (Number of Degrees Of Freedom), i.e. the number of vessel system’s key dimensions being virtually unlimited, and all vessel system’s critical components being coalesced/merged into 1D systems synthetic vector $\mathit{R}\left(\mathit{t}\right)$. Cargo vessel’s panel stress PSD (Power Spectral Density) within most stationary 3 h ocean/sea states containing two peaks (i.e. bimodal spectrum). WF (Wave Frequency) part was related to the wave-induced loads, and HF (High Frequency) part was caused by the transient load/resonant vibrations, known as springing/whipping. Subsequent research endeavours pertaining to the whipping phenomena need to scrutinize the whole dynamics of stress and accelerations in this study and discretely tackle the responses filtered at the low and high passes. The novel Gaidai multivariate risk assessment approach is applicable to higher dimensional cooperatively quasi-stationary, nonlinear dynamic systems. Key failure/damage mechanisms such as sudden shock damage, fatigue failures, corrosion-type cracks, or other ultimate collapse processes connected to the dynamic system of a cross-coupled cargo vessel may be introduced as additional dimensions or failure/damage modes to the dynamic system of interest if such measurements are available. The study’s long-term goal is to create a reliable code for ship route design that takes excessive loads and fatigue build-up into account. This is a crucial task, and future study may build on the suggested risk assessment strategy. It may be advised that further research to be done on other stochastic structural models for wind-wave environments in realistic oceanic regions. Potential research outcomes will greatly advance design practices of contemporary cargo vessels.

4. Discussion

The cargo ship on-board dataset was used to validate state-of-art multivariate risk assessment approach that this study offered. Most realistic choices for complex dynamic systems are out of reach, such as the direct MC simulation technique, lengthy laboratory tests, or, in the best-case scenario, on-board measurements. While the proposed Gaidai multivariate risk assessment methodology has no restrictions on the number of system dimensions, existing risk assessment methodologies utilized in naval design being mostly restricted to 2D systems. Traditional risk assessment techniques, pertaining to timeseries, lack the benefit of effectively managing systems with high dimensionality and cross-correlation across distinct critical components inside vessel hull’s dynamic system. The capability of the presented approach to examine the risk assessment of high-dimensional dynamic systems is a key benefit. Previous validation of the approach described in this study was done using a variety of simulation models, but only for one-dimensional system responses. Overall, highly accurate predictions have been made. This work aims to design a multivariate risk assessment approach that is resilient, easy to use, and general-purpose.

Regarding on-board measurements, confidentiality concerns mean that these datasets being rarely made accessible for public research presentations. Current study was able to make use of a genuine, sizable, raw, measured on-board dataset. The goal of this research is to support the creation of new DNV safety guidelines and standards for the modern marine transportation sector, which may have an impact on future cargo transport insurance and cargo vessel operating policies.

5. Conclusions

This study examined ship dynamic response timeseries and a synthetic windspeed data set that were collected on board during more than 70 transatlantic trips by 4400 TEU Panamax container ships in the 2010s. A detailed explanation of the theoretical basis for the suggested approach is provided. It should be noted that it is appealing to analyse the dynamic system risk assessment using direct measurement or large-scale MC simulations. The enormous dimensionality and complexity of dynamic systems, however, need the development of novel, reliable, yet accurate algorithms that can work with a small amount of accessible data while making the most use of it. Using both synthetic and real on-board measured ship panel stress datasets, the innovative technique was verified against the 4-parameter bivariate Weibull method bivariate method. Finally, a variety of engineering domains can employ the proposed technique. The example of naval architecture given does not restrict the domains in which Gaidai multivariate risk assessment approach methodology may be utilized.

For future studies, larger container vessels along with increased number of sensor readings at additional crucial points on deck panels can be studied. Larger TEU vessel hull whipping responses – which are concurrently captured at many vessel hull locations – to receive special attention. When discussing constraints and limitations of Gaidai multivariate risk assessment methodology, non-stationary dynamic vessel systems with underlying trends should be brought up. In the latter scenario, the Gaidai multivariate risk assessment approach can still be used once the underlying trend has been first determined and then removed from the raw information.

List of abbreviations

TEU	=	twenty feet equivalent unit
MC	=	Monte Carlo
AOG	=	Active Operator Guidance
CI	=	Confidence Interval
P	=	port, S – starboard side
DL	=	deck longitudinal
DLP	=	deck sensor on port side
AFV	=	vessel’s forward acceleration

Author contributions

All authors have contributed equally.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Data will be available on request from the corresponding author.

Additional information

Funding

The author(s) reported there is no funding associated with the work featured in this article.