Reliability-based assessment of ship hull girder ultimate strength

List of symbols

A p

net sectional area, in cm², of attached plating

A pE

effective area of the plate, in cm²

A s

net sectional area, in cm², of stiffener

a

ordinary stiffener spacing, in m

B

breadth of the ship, in m

B 0

maximum amplitude of initial imperfection for column and torsional buckling mode

b

distance of the transverse stiffener, in m

b E

effective width of the stiffener, in mm

C 1

correction factor due to randomness in yield stress

C 2

correction factor for yield stress reduction due to initial imperfection

C b

block coefficient of the ship

F M Wv

cumulative distribution of peak value of WIBM

F X sw

cumulative probability function of SWBM in one voyage

F X Wv

cumulative probability function of WIBM in the given period

f X sw

probability density function of SWBM in one voyage

h w

net web height, in mm

h we

effective height of the stiffener, in mm

L

length of the ship, in m

L FEM

length of the FEM model of stiffened panel

M SW

bending moment in still water condition (SWBM)

M SW,c

the characteristic value of SWBM

M t

total vertical bending moment

M Wv

wave-induced bending moment (WIBM)

M Wv,c

the characteristic value of WIBM

N

number of wave cycles encountered by the ship in the given period

n

number of voyages in a particular period

R eHp

minimum yield stress, in N/mm², of the plate

R eHs

minimum yield stress, in N/mm², of the stiffener

t f

thickness of the flange, in mm

t p

thickness of base plate, in mm

t w

thickness of the web, in mm

w 0

maximum amplitude of initial imperfection for local buckling mode

X SW,m

model uncertainty in SWBM estimation

X SW , st

statistical uncertainty in SWBM estimation

X Wv,envir

environmental uncertainty in WIBM estimation

X Wv,m

model uncertainty in WIBM estimation

X Wv,nl

uncertainty due to nonlinearity in WIBM estimation

X Wv,st

statistical uncertainty in WIBM estimation

List of Greek symbols

β

reliability index

β E

effective plate slenderness

γ

shape parameter of the extreme value distribution of SWBM

δ

scale parameter of the extreme value distribution of SWBM

ζ

scale parameter of the Weibull distribution of WIBM in each cycle

η

shape parameter of the Weibull distribution of WIBM in each cycle

μ X SW,st

mean value of statistical uncertainty of SWBM

μ X Wv,st

mean value of statistical uncertainty of WIBM

σ X SW,st

standard deviation of statistical uncertainty of SWBM

σ X Wv,st

standard deviation of statistical uncertainty of WIBM

σ C 1 , 2 , 4

critical stress in N/mm² for each respective buckling mode

σ CP

buckling stress of attached plate in N/mm²

σ CR 1

ultimate strength of stiffened panel under column buckling

σ CR 2

ultimate strength of stiffened panel under torsional buckling

σ CR 3

ultimate strength of stiffened panel made of flanged profiles under web local buckling

σ CR 4

ultimate strength of stiffened panel made of flat bars under web local buckling

σ CR 5

ultimate strength of plate buckling

Φ

edge function

ψ W

load combination factor

List of abbreviations

CSR

Common Structural Rules

HGUS

hull girder ultimate strength

IACS

International Association of Classification Societies

ISSC

International Ship and Offshore Structures Congress

SWBM

still water bending moment

WIBM

wave-induced bending moment

2.1 Reference model

The bulk carrier used for this present study was the ISSC-2000 Capesize bulk carrier [27]. The cross-section of its hull girder is shown in Figure 1. This type of bulk carrier has longitudinal strengthening except for the side shell between the top side tank and the bilge hopper tank, which is transversely stiffened [1].

Figure 1

Hull girder cross-section of the ISSC-2000 bulk carrier.

2.2 Case configuration

Mainly, the loading conditions of the bulk carriers are homogenous hold loading, AHL conditions, and cargo in two or more adjoining holds [21]. Among those loading conditions, the alternate loading condition is the most critical for the safe design of bulk carriers [1]. The alternate loading condition is applied for carrying high-density cargos in odd-numbered cargo holds, as shown in Figure 2. When the cargo is distributed in odd-numbered cargo holds, it gives double the cargo weight per hold of homogenous loading, precipitating the maximum shearing force of the hull girder caused by the opposite-acting vertical force of the hull/cargo weight vs buoyancy.

Figure 2

Cargo placement and shearing force direction in AHL condition.

The bulk carrier was assumed to be operated by crossing the North Atlantic and spending 23.5 days finishing 1 voyage [1]. Sea areas in the North Atlantic given by global wave statistics were used to evaluate the WIBM during operation [22], covering world sea areas 8, 9, 11, 15, 16, and 17, as shown in Figure 3 [34].

Figure 3

North Atlantic wave scatter data.

When conducting time-dependent reliability analysis, the case configuration was built from two parameters. The first was the duration of the operation, and the second was the duration of the ship spending in AHL condition. The ship operation duration was divided into short-term (1 voyage and 1 year) analysis and long-term (5–25 years with 5 years increments) analysis. The AHL condition was focused on two assumptions. First, the ship was assumed to be in AHL condition for 50% of voyages, and second, the ship was assumed to be in AHL condition for 20% of voyages. The first assumption was that the ship only had one port for unloading at its final destination. The second was to cover the possibility of the ship having two ports for unloading. From the loading port to the first unloading port, the ship was in homogenous loading condition. From the first unloading port to the last unloading port, indicating 20% of the voyage, the ship was in AHL condition. The duration of the ship in the AHL condition affects the number of waves encountered by the ship during operation and then influences the statistical uncertainty of the WIBM. The case configuration is shown in Table 1 and analyzed in both sagging and hogging conditions.

3.1 Incremental-iterative method (Smith’s method)

In Smith’s method, the ultimate strength of the hull girder is calculated at the transverse sections between adjacent transverse webs. A hull girder cross-section is first divided into a set of independent elements. In each incremental step, the incremental curvature is imposed, and the bending moment is generated. Curvature increments correspond to rotation angles of the transverse section of hull girders around their horizontal neutral axis. Then, the corresponding strain and stress are calculated for each element. This procedure is repeated for every step with an updated neutral axis by establishing force equilibrium over the whole transverse section. The corresponding moment is calculated by summing the contributions of all elements. Finally, the maximum HGUS is obtained from the peak value of the curve plotting the vertical bending moment vs the curvature [28]. In the hogging condition, the structural elements above the neutral axis are lengthened while the elements below the neutral axis are shortened, and vice-versa in the sagging condition.

In order to calculate the HGUS, the hull girder cross-section is divided into particular members, where each member contributes to the hull girder strength. The structural members are categorized into stiffener elements, stiffened plate elements, and hard corner elements. The details on this structural element classification and the formulations describing the load-shortening curve for each failure mode of the structural elements are referred to in [21] as listed in Eqs. (1)–(5).

Eq. (1) calculates the stiffener element’s load shortening curve in the column buckling mode, Eq. (2) in the torsional buckling mode, and Eqs. (3) and (4) in the web local buckling mode for flanged profiles and flat bars. These three load-shortening curves are computed for each stiffener element. The maximum stress values for the various failure modes are compared, and the ultimate strength of the stiffener element is determined as the minimum value. In the case of a stiffened plate element, the plate buckling mode of Eq. (5) is used. More in-depth explanations can be found in the IACS-CSR code [21].

Considering its stochastic behavior, the yield stress input data for R eHs and R eHp in Eqs. (1)–(5) can be described as follows:

C 1 is the correction factor due to randomness in yield stress, and C 2 is the correction factor due to randomness in initial imperfection.

3.2 Uncertainties in HGUS

3.2.2 Initial geometric imperfection

The proposed scatter data of initial geometric imperfection in ref. [19] was adapted as shown in Figure 4.

Figure 4

Analyzed scatter data of initial imperfection.

The maximum amplitude of the initial geometric imperfection was calculated as follows:

(7) w 0   = 0.1   β p 2 t ( mm ) for β p 2 t < 60 ,

  w 0   = 6   ( mm )  for  β p 2 t > 60 .

The maximum amplitude of initial geometric imperfection was defined as a function of the slenderness of plating, calculated as 10 3 s t p R e H E . The mean value was set as half of the maximum amplitude of the initial geometric imperfection, and the standard deviation was set to be half of the mean value. The value of β ² t for all stiffener members on the hull girder cross-section was calculated, resulting in all stiffener members being located inside the blue line of Figure 4. Thus, to strengthen the focus of the analysis, the region inside the blue line was particularly investigated. The increment of the geometric imperfection was gradually increased from a small geometric imperfection to 100% of the calculated maximum geometric imperfection.

As stated in Eq. (6), the effect of initial geometric imperfection on the ultimate strength of a hull girder was imposed in terms of yield strength reduction of the respective elements. For this purpose, finite element method (FEM) was first employed to determine the effect of initial geometric imperfection on stiffener elements. However, using FEM to analyze all sizes of the stiffeners attached to the hull girder is time-consuming. It is necessary to obtain a simplified formula for forecasting the decrease in elements’ yield strength due to geometric imperfection for various values of element inertia.

In this calculation stage, it is convenient to introduce idealized models for representing the rest of the stiffener elements. Calculation accuracy depends heavily on how representative the calculation models are of all stiffener elements in the cross-section of the hull girder. Concerning this matter, the outer borders of the area inside the blue line in Figure 4 were set to be the idealized models. Furthermore, the stiffener with a value of β ² t near 60 was also chosen as another idealized model. Thus, three idealized models were selected, namely the light model, the medium model, and the heavy model, respectively. The scantlings of these idealized models are shown in Figure 5.

Figure 5

Scantlings of idealized models.

The effect of initial geometric imperfection on the ultimate strength of the stiffened panel is influenced by not only the maximum amplitude of initial imperfection but also the buckling modes of failure, including the column buckling mode, the lateral-torsional buckling mode, and the web buckling mode [21,22]. Thus, the initial geometric imperfection should be imposed concerning all possible buckling modes. Following this, the total initial geometric imperfection was defined as the sum of all the imperfections of the buckling modes, as shown in Figure 6.

Figure 6

Imperfection modes: (a) column buckling mode; (b) torsional buckling mode; and (c) local buckling mode.

The formulations for generating the initial geometric imperfection are listed in Table 4.   B 0 is the maximum initial geometric imperfection for the column and lateral-torsional buckling as B 0 = L FEM / 1 , 000 . C 0 is the maximum initial geometric imperfection for local buckling mode, calculated as w 0 in Eq. (7), which depends on the slenderness of plating β p [28,29,30].

Table 4

Deformation of initial geometric imperfection

Imperfection mode	Formulation
Column buckling	w 0 column = B 0 sin π x i a
Torsional buckling	w 0 torsional = B 0 z i h s sin π x i a
Local buckling	w 0 local = C 0 sin m π x i a sin π y i b

To ensure the accuracy of the ultimate strength calculation of stiffened panels using FEM, the boundary conditions and the extent of the model must be carefully considered. The model proposed by Yao et al. [28], also used in [29], was adopted. The model is referred to as the 1/2 + 1 + 1/2 bay and 1/2 + 1 + 1/2 span model based on the covered bay and span (Figure 7).

Figure 7

1/2 + 1 + 1/2 bay and 1/2 + 1 + 1/2 span model.

For the boundary condition, the periodic continuous condition was considered for both longitudinal and transverse edges by setting the displacements along all edges uniform, as shown in Figure 8. To simplify the calculation, the transverse frames were not modeled. Instead, their role was considered as a deflection constraint applying to the plate and the stiffener along the line of the transverse frames [29,30].

Figure 8

Boundary condition of FEM.

3.3 HGUS calculation results

3.3.1 Results of the initial geometric imperfection

This section aims to obtain the relationship between the inertia of the area of stiffened panel and its normalized ultimate strength for various initial imperfection conditions. The normalized ultimate strength refers to the comparison value between the ultimate strength obtained by FEA and IACS-CSR. For the first step, the load-shortening curve for all idealized models was calculated using Eqs. (1)–(4), resulting in Figures 9–11.

Figure 9

The load-shortening curve of the light model based on IACS-CSR.

Figure 10

The load-shortening curve of the medium model based on IACS-CSR.

Figure 11

The load-shortening curve of the heavy model based on IACS-CSR.

The calculation results for all buckling modes were compared, and the ultimate strength of stiffened panels was determined as its minimum value. From Figures 9–11, it is shown that all idealized models collapsed in column buckling mode. After obtaining the ultimate strength of a stiffened panel using IACS-CSR, the ultimate strength for various initial imperfection conditions was calculated using the finite element method. The results of FEM are shown in Figure 12.

Figure 12

Ultimate strength for various conditions of initial imperfection by FEM.

δ is the maximum calculated initial geometric imperfection ( B 0 and C 0 in Table 4) and δ ₀ is the maximum initial geometric imperfection for particular cases. From Figure 12, it can be seen that the value proposed by IACS-CSR is close to the finite element results for the initial imperfection (δ ₀/δ) of 50% in the light model and 75% in the medium and heavy models.

The next step was normalizing the FEM results in Figure 12 with the IACS-CSR results (σ/σ _CSR). The normalized value was then used as a correction factor C 2 in Eq. (6). The correction factor C 2 for the idealized models in various conditions of initial imperfection is plotted in Figure 13 as a straight line, and its regression equation for each increment of initial imperfection is generated accordingly. The correlation between the correction factor and the inertia of the stiffener element predicted by the regression equation is plotted in the dashed line. Comparing the straight and dashed lines, it can be seen that the regression equation accurately predicts the correction factor for each respective initial imperfection.

Figure 13

Correction factor due to initial geometric imperfection.

3.3.2 HGUS

Considering the uncertainties in material properties, geometric properties, initial geometric imperfection, and corrosion behavior and applying Monte Carlo simulation for 500 simulations per case, the distribution of HGUS versus the operation years was plotted in Figure 14 for hogging and in Figure 15 for sagging conditions, respectively. Roughly, HGUS decreases linearly in the case of a standard corrosion rate but follows a polynomial trend line in the case of a severe corrosion rate. In severe conditions, after 15 years of operation, some stiffeners had already exceeded the corrosion allowance. As stiffeners were being replaced, the decrease in HGUS for the following years slowed down.

Figure 14

HGUS vs operation years in hogging condition.

Figure 15

HGUS vs operation years in sagging condition.

4.1 Statistical description and load combination

The uncertainties associated with SWBM are caused by various factors, including loading patterns, trade routes, cargo density, and human factors. In conventional ships, SWBM generally follows a normal distribution [1]. The SWBM is defined as follows:

X SW,st is statistical uncertainty, X SW,m is model uncertainty, and M SW,c is the characteristic value of SWBM based on IACS-CSR [21] as follows:

where

In the assessment of WIBM, the considered uncertainties are statistical uncertainties, model uncertainties, environmental uncertainties, and uncertainty due to the nonlinearity effect, which can be mathematically described as follows:

X Wv,st , X Wv , m , X Wv,nl , and X Wv,envir are the statistical, model, nonlinearity, and environmental uncertainties, and X Wv , c is the characteristic value of the WIBM. Corresponding to an exceedance probability of 10⁻⁸, the characteristic value of wave-induced hogging and sagging bending moments are assessed in kN m as follows:

where the peaks in each wave cycle distribute according to the Weibull distribution:

where ζ and η are the scale and shape parameters. The shape parameter typically ranges from 0.9 to 1.1 for bulk carriers [33]. Here, by defining the shape factor as 1, the Weibull distribution was assumed to reduce to an exponential distribution. The scale parameter can be expressed in relation to the characteristic value as follows:

Because SWBM and WIBM are random, their maximum values rarely coincide. As a result, the total vertical bending moment is less than the sum of SWBM and WIBM. To deal with this phenomenon, the total bending moment was calculated as Eq. (17), with the load combination factor ψ W introduced as 0.97 [1]:

4.2 Uncertainties in SWBM

As Eq. (8), there are two uncertainties in SWBM calculation: model and statistical uncertainties. X S W , m is the model uncertainty in SWBM prediction caused by varieties of loading patterns and cargo densities. It follows a normal distribution with a mean value of 1 and a standard deviation of 0.1.

X S W , s t denotes the statistical uncertainty related to the number of voyages in the reference period. It follows the Gumbel distribution, with the mean value and the standard deviation calculated as follows:

where γ is the shape parameter and δ is the scale parameter of the extreme value distribution of the SWBM. It can be approximately assumed to follow a Gumbel distribution calculated from the initial normal distribution of SWBM in 1 voyage [36] as follows:

where n is the number of voyages in a particular period, X sw is a random variable denoting SWBM in the AHL condition, X sw follows the normal distribution with a mean value of 0.7, a standard deviation of 0.1 in 1 voyage duration [1,36], F X sw − 1 is the inverse cumulative probability function of SWBM in 1 voyage, and f X sw is the probability density function of SWBM in 1 voyage.

4.3 Uncertainties in WIBM

The first is the statistical uncertainty X Wv,st . Assuming that the distribution of X Wv,st converges to the Gumbel extreme value distribution, the cumulative distribution function of the maximum WIBM in N wave cycles can be obtained as follows:

where λ and α are parameters of the Gumbel distribution, which can be related to the initial Weibull distribution of Eqs. (15) and (16) as follows:

where N is the number of wave cycles encountered by the ship in the given period. The mean value and standard deviation of the Gumbel distribution can be obtained as follows:

The second is the model uncertainty X wv,m caused by the use of linear theory in the long-term prediction of WIBM. It can be assumed to follow the normal distribution with a mean value of 1.0 and a standard deviation of 0.1 [1].

Environmental uncertainty X wv,envir is basically a bias factor for considering variations in wave scatter diagrams. It is a ratio between WIBM based on other wave scatter data and WIBM based on IACS data. Since this study uses the IACS wave data as the reference and the IACS-CSR to calculate the characteristic design value, the bias factor is set as 1.0.

The last is uncertainty due to the nonlinearity effect X wv,nl . For ships with large block coefficients, such as bulk carriers, the nonlinear hydrodynamic effects cannot be denied. Since the characteristic WIBM by IACS-CSR used here accounts for the effect of nonlinearity, it can be considered normally distributed with a mean value of 1.0 and a standard deviation of 0.1 [1].

4.4 Calculation result

In order to analyze how the total loads change over the operation period, the relationship between the increase in loads and the years of operation was plotted (see Figure 16 for the hogging condition and Figure 17 for the sagging condition).

Figure 16

Total loads vs operation year in hogging condition.

Figure 17

Total loads vs operation year in sagging condition.

The dashed line indicates the mean value of the total loads in the hogging and sagging conditions for an AHL duration of 50% and the straight line for an AHL duration of 20%. The lines marked with a triangle and a square show the range of majority occurrences. From Figures 16 and 17, the difference in the mean value of the total loads between the AHL duration of 20% and the AHL duration of 50% is relatively constant for all years of operation in both hogging and sagging conditions.