Ultimate Limit State Function and Its Fitting Method of Damaged Ship under Combined Loads

Abstract

The ultimate limit state function is one of the premises for the assessment of structure strength and the safety of ships under severe conditions. In order to study the residual strength of damaged ships under the combined load of vertical and horizontal bending moments acting on the hull girder, the ultimate limit state function of a damaged ship under combined load, and its fitting methods are investigated in this paper. An improved Smith Method is adopted to calculate the residual load carrying capacity of damage ships, where the rotation and translation of the neutral axis of the damaged cross-section are obtained using a particle swarm optimisation method. Because the distribution curve of the residual load carrying capacity of a damaged ship under combined load is asymmetric, the application of traditional explicit polynomial fitting methods results in poor accuracy. In this study, a piecewise weighted least square fitting method is adopted so as to guarantee the continuity in the transitions, and a method is proposed for fitting the ultimate limit state function of a damaged ship under combined load. Calculations of the residual strength show that the improved Smith Method is more accurate than the original Smith Method for the accurate position of the neutral axis. The error analysis of the fitting methods shows that the ultimate limit state function that is fitted using a piecewise weight least square method is more accurate.

1. Introduction

Ships may encounter various types of dangers at sea; for instance, list due to stranding or hull damage due to collision. The Ro-Ro ship Modern Express underwent an accident in 2016, and floated at sea for more than one week with a list angle of 51 degrees. Compared to ships floating in the upright position, the hydrostatic load acting on a damaged ship is very different. When a ship is in a listed position, the horizontal bending moment, in addition to the vertical bending moment, must be also considered in the body-fixed coordinate system. If the ship hull is also damaged, the strength of the hull girder will be impaired, and the ship will be in a more dangerous situation. Researches have shown that the horizontal load acting on a damaged ship may be as large as 1.73 times the vertical load [1,2]. For this reason, the assessment of ship strength is more than just about dealing with the vertical bending moment; in fact, the ultimate limit state function of the ship under the action of combined loads must be assessed.

To investigate the ultimate limit state function, the ultimate strength must be calculated first. The existing methods for calculating the ultimate strength of the hull girder can be categorised into two groups: analytical methods and progressive collapse methods [3]. The analytical methods include the elastic analysis method (such as the initial yielding bending moment method) and fully plastic bending moment method (Caldwell’s method [4]). The practices show that these methods are useful for predicting the ultimate strength and residual strength of the hull girder in the early design stage. Based on investigations of the collapse procedure of a ship undergoing the bending moment of a hull girder, some progressive collapse analysis methods have been proposed, such as nonlinear finite element analysis (including the Idealised Structural Unit Method) [5,6] and Smith’s method [7]. The nonlinear finite element analysis, which can simultaneously consider yielding, buckling and the mixed failure behaviour of structural elements constituting the hull girder, is recognised as the most reliable approach for calculating the ultimate/residual strength. However, its application is limited due to the huge modelling and computing time, and also the requirement of the engineer’s experience [8,9]. Meanwhile, a progressive collapse method was proposed by Smith [7], by which, the cross-section of concern is discretised into units of stiffened panels, then the curvature of the hull girder is increased step by step, and the nonlinear behaviour of each unit is obtained by analysing its stress–strain relationship so as to obtain the section bending moment at different curvatures of the hull girder. This method has been widely used in the community, and is referred to as the Smith Method [10]. Hence, the simplified progressive collapse method is adopted for the ultimate strength analysis in this study.

The accuracy of the Smith Method relies on the stress–strain relationship and the calculation of the neutral axis. There exist quite a few researches on the stress–strain relationship. Hughes et al. [11] investigated the failure behaviour of various types of stiffened panels using an incremental iteration method, and the comparison with experiment showed fair accuracy. Yao [12] derived the stress–strain relationship for the coupled flexural–torsional behaviour of angle bar stiffeners. HCSR [10] also provided various stress–strain relationships for the calculation of the ultimate strength of hull girders. However, the relationships provided by HCSR were derived using the compression and tension equilibrium in the structure above and below the neutral axis that is parallel to the water line, and they cannot be applied to asymmetric cross-sections due to inclination and/or hull damage. In the case of an asymmetric stress–strain, the neutral axis does not only translate, but also rotate. Fujikubo et al. [13] found that the influence of the rotation of the neutral axis on the accuracy of residual strength may be as large as 8%. In order to account for the rotation, Joonmo et al. [14] proposed a force vector equilibrium method where both the force and moment at the cross-section are in equilibrium.

Although theoretically, the instantaneous position of the neutral axis can be accurately determined by the equilibrium both in force and force vector, a proper numerical method is required for the convergence of a solution. Li et al. [15] proposed a linear search method to obtain the solution. However, the translation and rotation are dealt with separately in the iteration process, and certain experience is needed to obtain reasonable results.

The solution to the instantaneous position of the neutral axis is a problem of seeking the optimum solution, and for which there exist a number of algorithms, such as Particle Swarm Optimisation (PSO) [16], ant colony optimisation [17], genetic algorithm [18] and simulated annealing [19]. Among these algorithms, the PSO algorithm was inspired by the social behaviour of bird swarms. Massless and volumeless particles are adopted to represent individuals of the swarm, and each individual is a solution. The motion of each of these particles obeys a simple law, and the global optimum solution is searched in the solution space, accounting for the influences between the particles. This algorithm is widely applied in many engineering areas for its simple implementation, fast convergence, and excellent accuracy. Li et al. [20] applied the PSO algorithm to trace the instantaneous position of the neutral axis, and both the translation and rotation are dealt at the same time. Comparison of results showed that the accuracy is satisfactory.

Among the researches on the fitting of ultimate limit state function under multiple loads, mostly the combined load of vertical and horizontal bending moments [21,22,23], the following fitting equation is often used.

$${\gamma }_1{\left(\frac{M_V}{M_{UV}}\right)}^{\alpha }+{\gamma }_2{\left(\frac{M_H}{M_{UH}}\right)}^{\beta }=1,$$

(1)

where ${\gamma }_1$, ${\gamma }_2$, $\alpha$ and $\beta$ are the fitting parameters of the function; $M_V$ and $M_H$ the vertical and horizontal bending moment of the ship in the ultimate state, respectively; $M_{UV}$ the ultimate vertical bending moment when the ship is in a sagging or hogging state; and $M_{UH}$ the ultimate horizontal bending moment when included to port or starboard. These parameters vary from one to another, and

Table 1. Fitting parameters used by researchers for intact ships.

Contributor		${\mathit{\gamma }}_1$	${\mathit{\gamma }}_2$	$\mathit{\alpha }$	$\mathit{\beta }$
Gordo et al. [21]		1	1	$1.5<\alpha =\beta <1.66$
Mansour et al. [22]	If $\left|\frac{M_V}{M_{UV}}\right|>\left|\frac{M_H}{M_{UH}}\right|$	1	0.8	1	2
	If $\left|\frac{M_V}{M_{UV}}\right|<\left|\frac{M_H}{M_{UH}}\right|$	0.8	1	2	1
Khan and Das [23]		1	1	Not fixed	Not fixed

shows some of them used by researchers for intact ships.

It can be seen in Table 1 that attempts have been made to obtain a function to describe the coupling relationship between the two bending moments. The fitting is satisfactory when the samples of bending moments bear a negative relationship with the absolute values; however, if the samples bear a positive relationship, it will be very difficult to obtain a fitting function with desired accuracy, as shown in Figure 1, where the trend of the samples starts to change when approaching the $x$ axis, and increasing the number of samples is of little help, as shown in Figure 1.

In order to improve the fitting accuracy of ultimate limit state function, Shahid [24] proposed a polynomial fitting method for the case of damaged ships. The polynomial fitting method achieves a quadratic fitting function using multivariable regression. However, it has been shown that there are still some errors in the results obtained with this method. Thus, a new fitting method is desired to accurately describe the relationships between the vertical and the bending moment.

In this study, in order to derive a simple but accurate fitting method, the accuracy of the polynomial fitting method is improved by the means of a piecewise weighted least squared method. With this method, the sample points are divided into groups and fitted separately, then these functions are adjusted by ensuring the continuity at the boundaries, and finally a set of explicit fitting functions are obtained.

2. PSO-Based Smith Method for Ultimate Bearing Capacity Analysis of Damaged Hull Girder

2.1. Smith Method and the Neutral Axis of Asymmetric Section of Damaged Ship

The Smith Method is widely used for the assessment of the ultimate strength of ships, by which, the cross-section of ship hull girder is discretised into units of stiffened panels and hard corner, and the ultimate strength is obtained by an increase of the curvature of the hull girder step by step. The solving process of the Smith Method is introduced in HCSR [10] in detail. The accuracy of the Smith Method relies on the stress–strain relationship and the calculation of the neutral axis.

In most of the cases of an intact ship in an upright floating position, there exist three symmetries: symmetry in geometry, symmetry in material, symmetry in load. Intact ships are usually symmetrical in geometry and material, and the most common two types of asymmetry are:

(1)

Geometrical asymmetry of damaged ships;

(2)

Load asymmetry due to an uneven distribution of goods, the flooding when the hull is damaged, or due to combination of the horizontal and vertical wave bending moments.

To apply the Smith Method, the position of the neutral axis must be determined. The convergence factor provided by HCSR for determining the position of the neutral axis takes the form

$$\frac{\left|F_c-F_t\right|}{\left|F_c+F_t\right|}\le {\xi }_T,$$

(2)

where $F_c$ is the sum of the compressive forces of all units, $F_t$ the sum of the tensile forces of all units, and ${\xi }_T$ is the convergence factor for force equilibrium, which is usually smaller than 0.01.

It can be seen that Equation (2) is applicable only when all the three symmetries are present. Absence of any of the three symmetries would result in an asymmetric cross-section. For instance, in the case of a damaged hull, the cross-section is asymmetric due to flooding or loss of goods, and the force acting on it will be as shown in Figure 2. The rotation of the neutral axis is not accounted for in the original HCSR rules, and application of it with any correction will result in a wrong instantaneous position of the neutral axis in this case, which will then introduce errors into the subsequent calculation of the stresses in the units, and eventually lead to an inaccurate solution to the ultimate bending moment.

A new equilibrium criterion is required to determine the position of neutral axis of the asymmetric cross-section. According to the research of Choung et al. [14], an additional force vector equilibrium criterion, in addition to the force equilibrium shown in Equation (2), is to be satisfied:

$$\frac{\pi }{2}-{\xi }_R\le {\cos }^{-1}\left(\frac{\stackrel{\rightharpoonup }{F}\cdot \stackrel{\rightharpoonup }{M}}{\begin{array}{cc}\left|\stackrel{\rightharpoonup }{F}\right|& \left|\stackrel{\rightharpoonup }{M}\right|\end{array}}\right)\le \frac{\pi }{2}+{\xi }_R,$$

(3)

where $\overrightarrow{F}$ is the sum of the forces acting on the cross-section, $\overrightarrow{M}$ the total section bending moment, and ${\xi }_R$ the force vector convergence factor (usually smaller than 0.1°).

2.2. The PSO-Based Smith Method

The instantaneous position of neutral axis can be obtained once the equilibrium criteria are met. Since both the transition and rotation need to be accounted for, an appropriate algorithm is required to obtain the solution with fair accuracy. This, actually, is a problem to search for the optimum solution in a given space [20].

In this study, the Particle Swarm Optimisation (PSO) method is used to obtain the solution. The PSO method approaches the optimum solution by collaborating and sharing information between the individuals of a group. Each particle has two properties: speed and location. Speed refers to the speed at which the particle is moving, and location indicates the direction of the motion. The optimum solution obtained for an individual particle is called individual best, while the optimum solution of all the particles is referred to as the global best. According to the individual best and global best, a new calculation is executed to update the speed and location of all the particles. This process is repeated until the convergence criteria are satisfied.

By the Smith Method, the position and rotation of the neutral axis will vary every time the curvature is increased, thus the two motions can be regarded as the two dimensions of the search space. Then the PSO can be applied to seek for the optimum global solution to the objective function. Since the translation and rotation vary at the same time instead of one by one, the actual instantaneous position of the neutral axis at a particular curvature can be obtained in a more direct manner. The flowchart of the PSO-based Smith Method, referred to as the improved Smith Method hereafter, is shown in Figure 3.

3. Fitting for the Ultimate Limit State Functions of Damaged Ships under Combined Bending Moments

3.1. Explicit Ultimate Limit State Functions

The curve of the ultimate limit state function indicates the ultimate vertical and horizontal load carrying capacity under combined load. It is difficult to accurately describe the curve in an explicit form. For this reason, the function is usually obtained by fitting with a number of samples.

Since it is a load combined from bending moments in two directions, i.e., vertical and horizontal, the rotation of the cross-section may be achieved by rotating the coordinate system. In this study, the improved Smith Method is applied to obtain the ultimate bending moment at each rotation. The corresponding ultimate rotation can then be decomposed into components in the axis directions to obtain the vertical and horizontal ultimate bending moment, as shown in Figure 4, where $Y_{0}O{Z}_0$ is the global coordinate system; $Y_{1}O{Z}_1$ the body-fixed coordinate system; $\theta$ the rotation; $M_0$ the total bending moment; and $M_V$, $M_H$ the vertical and horizontal bending moment respectively, which can be calculated from $M_0$

$$\{\begin{array}{l}{M}_V=M_0\cos \theta \\ M_H=M_0\sin \theta \end{array},$$

(4)

The double hull Very Large Crude Container (VLCC) model provided by ISSC [25] is used in this study. In order to guarantee the accuracy, the calculation is performed at every 10 degrees of rotation, and the result obtained for the total ultimate bending moment is decomposed into a vertical bending moment and a horizontal bending moment. The definitions of the horizontal bending moment components are shown in Figure 5. The vertical bending moment $y_i$ and the horizontal bending moment $x_i$ are obtained for each sample point.

The following equation [23,26] is often used for the explicit fitting.

$${\left(\frac{M_H}{M_{UH}}\right)}^{{\alpha }_1}+{\left(\frac{M_V}{M_{UV}}\right)}^{{\alpha }_2}=1,$$

(5)

where $M_{UV}$ and $M_{UH}$ are the ultimate bending moment resulted from pure vertical and horizontal bending in the global coordinate system, respectively. It can be written in the form

$$x^{{\alpha }_1}+y^{{\alpha }_2}=1,$$

(6)

where $x=M_H/M_{UH}$ and $y=M_V/M_{UV}$.

It can be seen that this is a nonlinear function, and which will be fitted using a nonlinear least square method in this study. This fitting method seeks the best match for the fitting function and determines the unknown coefficients via the least square of errors.

3.2. The Polynomial-Fitting-Based Ultimate Limit State Functions

A polynomial fitting method was proposed by Shahid [24] for the fitting of relationship between bending moments in different directions. Equation (7) shows the general form of the polynomial fitting function:

$$\overline{Y}=C_0+{\displaystyle \sum _{i=1}^n{C}_i{X}_i}+{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^n{C}_{ij}{X}_i}}{X}_j+\epsilon ,$$

(7)

where $\overline{Y}$ is the estimate of the fitting function for the $n$ random variables $X_i$, $C_0$ an independent regression coefficient, $C_i$ the coefficient of the first order term, $C_{ij}$ the coefficient of the cross-quadratic term, and $\epsilon$ the error between the estimate and the real response.

In this study, the third-order polynomial fitting function is adopted.

$$y=a+bx+c{x}^2+d{x}^3,$$

(8)

where $x=M_H/M_{UH}$ and $y=M_V/M_{UV}$. Since concave curves, i.e., multiple-value functions, may appear in some cases, the following form is used instead.

$$x=a+by+c{y}^2+d{y}^3,$$

(9)

Determination of the unknown coefficients is the key to the solution when applying the method, and, for which, the least square method is usually used.

3.3. The Improved Polynomial-Fitting-Based Ultimate Limit State Functions

For highly nonlinear systems, fitting with a single function usually results in poor accuracy, and increasing the order of the fitting function does not always help. Since the local variations of the samples are relatively small, a piecewise regression method is adopted to improve the accuracy of the polynomial fitting function, i.e., piecewise functions are used to describe the response relationship.

The accuracy can be improved by fitting with piecewise functions, but there still might be some large errors in some cases. For this reason, weights are introduced in order to obtain the optimum fit. To be specific, weights are used to improve the local fitting accuracy of each piece. The unknowns in Equation (8) can be calculated as follows.

$${\left\{N\right\}}^j={\left(X^{T}wX\right)}^{-1}{X}^{T}wY,$$

(10)

where

$$\left\{N\right\}={\left\{a,b,c,d\right\}}^T,$$

(11)

$$\left[X\right]=\left[\begin{array}{llll}1& X_1& {X_1}^2& {X_1}^3\\ 1& X_2& {X_2}^2& {X_2}^3\\ ⋮& ⋮& ⋮& ⋮\\ 1& X_m& {X_m}^2& {X_m}^3\end{array}\right],$$

(12)

$$\left\{Y\right\}={\left\{y_1,y_2,\dots ,y_m\right\}}^T,$$

(13)

where $w$ is the weight.

Then the weights are modified repeatedly until the accuracy requirement is satisfied. Once the fitting is done for all of the pieces, the values of the fitting functions at both ends of the piece, $y_{i1}$ and $y_{i2}$, are obtained, and so are the slopes $y_{i1}^{\prime }$ and $y_{i2}^{\prime }$. In order to ensure that the two fitting functions of the two adjacent pieces have the same curvature at the joint, the following boundary conditions are applied:

$$y_{i2}=y_{(i+1)1}=y_i,$$

(14)

$$\{\begin{array}{l}{y}_{i1}^{\prime }=(y_{(i-1)2}^{\prime }+y_{i1}^{\prime })/2\\ y_{i2}^{\prime }=(y_{i2}^{\prime }+y_{(i+1)1}^{\prime })/2\end{array},$$

(15)

Because it is an enclosed curve, the following boundary conditions are also applied

$$\{\begin{array}{l}{y}_{(n+1)1}=y_1\\ y_{(n+1)1}^{\prime }=y_{11}^{\prime }\end{array}.$$

(16)

The fitting function of each piece can be recalculated once the function values and the slopes at the ends of each piece are determined

$$\{\begin{array}{l}{y}_1=a+b{x}_1+c{x}_1^2+d{x}_1^3\\ y_1^{\prime }=b+2c{x}_1+3d{x}_1^2\\ y_2=a+b{x}_2+c{x}_2^2+d{x}_2^3\\ y_2^{\prime }=b+2c{x}_2+3d{x}_2^2\end{array},$$

(17)

4. Numerical Simulations and Results Analysis

4.1. The Details of VLCC

The cross-section of a double hull VLCC provided by ISSC [25] is adopted for the investigation; see Figure 6. The types and dimensions of the stiffeners are listed in

Table 2. The stiffeners and dimensions of the double hull VLCC [25].

Stiffener ID	Dimensions (mm)	Type	Yield Stress (MPa)	Stiffener ID	Dimensions (mm)	Type	Yield Stress (MPa)
1	300 × 90 × 13/17 IA	Angle Bar	313.6	25	250 × 90 × 12/16 IA	Angle bar	313.6
2	350 × 100 × 12/17 IA	Angle Bar	313.6	26	450 × 11 + 150 × 22	T-bar	352.8
3	400 × 100 × 11.5/17 IA	Angle bar	313.6	27	450 × 11 + 150 × 19	T-bar	352.8
4	400 × 11 + 150 × 12	T-bar	313.6	28	450 × 11 + 150 × 16	T-bar	352.8
5	400 × 11 + 150 × 14	T-bar	313.6	29	450 × 11 + 150 × 14	T-bar	352.8
6	450 × 11 + 150 × 12	T-bar	313.6	30	450 × 11 + 150 × 12	T-bar	352.8
7	400 × 11 + 150 × 14	T-bar	313.6	31	450 × 11 + 150 × 14	T-bar	352.8
8	450 × 11 + 150 × 16	T-bar	313.6	32	400 × 100 × 11.5/16 IA	Angle bar	352.8
9	450 × 11 + 150 × 19	T-bar	313.6	33	350 × 100 × 12/17 IA	Angle bar	352.8
10	450 × 11 + 150 × 22	T-bar	313.6	34	300 × 90 × 13/17 IA	Angle bar	352.8
11	450 × 11 + 150 × 25	T-bar	313.6	35	850 × 17 + 150 × 19	Angle bar	352.8
12	500 × 11 + 150 × 28	T-bar	313.6	36	250 × 90 × 12/16 IA	Angle bar	352.8
13	500 × 11 + 150 × 30	T-bar	313.6	37	300 × 90 × 12/16 IA	Angle bar	352.8
14	500 × 11 + 150 × 32	T-bar	313.6	38	400 × 11 + 150 × 14	T-bar	352.8
15	500 × 11 + 150 × 34	T-bar	313.6	39	450 × 11 + 150 × 12	T-bar	352.8
16	550 × 12 + 150 × 30	T-bar	313.6	40	450 × 11 + 150 × 14	T-bar	352.8
17	550 × 12 + 150 × 25	T-bar	313.6	41	450 × 11 + 150 × 16	T-bar	352.8
18	350 × 100 × 12/17 IA	Angle bar	313.6	42	450 × 11 + 150 × 19	T-bar	352.8
19	550 × 12.5 + 150 × 32	T-bar	352.8	43	450 × 11 + 150 × 22	T-bar	352.8
20	500 × 11.5 + 150 × 30	T-bar	352.8	44	450 × 11 + 150 × 25	T-bar	352.8
21	500 × 11.5 + 150 × 28	T-bar	352.8	45	450 × 11 + 150 × 28	T-bar	352.8
22	500 × 11 + 150 × 25	T-bar	352.8	46	500 × 11 + 150 × 25	T-bar	352.8
23	450 × 11 + 150 × 28	T-bar	352.8	47	500 × 11 + 150 × 28	T-bar	352.8
24	250 × 12.5	Flat bar	313.6	48	230 × 12.5	Flat bar	313.6

, and the Young’s modulus and the Poisson ratio are 206 GPa and 0.3, respectively. Typical cases of hull damage due to collision and stranding for the simulation are determined according to the HCSR, and the definitions of the damage are shown in

Table 3. Location and dimensions of the damage [10].

Dimensions(m)	Damage Due to Stranding	Damage Due to Collision
Height, D	Min(B/20, 2)	0.60D
Width, B	0.60B	B/16

and Figure 7.

4.2. Validation of PSO-Based Smith Method

In order to validate the PSO-based Smith Method (Method 3), the ultimate strengths of the three cases shown in

Table 4. Cases for validating the PSO-based Smith Method.

Case No.	Hull Intactness	State of Floating
Case 1	Intact	upright
Case 2	Damage due to collision	Inclined to Portside by ${30}^{\circ }$
Case 3	Damage due to stranding	Inclined to Portside by ${30}^{\circ }$

are calculated, and the results are compared to those obtained with the original Smith Method (Method 1) and the linear search-based Smith Method (Method 2).

The ultimate strengths of the double hull VLCC in sagging and hogging conditions of each case are calculated with the three Smith methods. The results for the ultimate strengths are shown in Figure 8, Figure 9 and Figure 10, and those for the positions of the neutral axis are shown in Figure 11, Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16. The ultimate sectional bending moments in the sagging and hogging condition of CASE 1 are compared with literature in

Table 5. Comparison of ultimate strength obtained by different contributors.

Source	Contributor	Methodology	Ultimate Strength (MN·m)
			Hogging	Sagging
ISSC 2000 Report [25]	Chen	ISUM	27.40	24.33
	Cho	Smith Method	28.66	20.80
	Yao	Smith Method	28.88	20.42
	Rigo (1)	Smith Method	28.31	19.57
	Rigo (2)	Modified P-M	25.61	24.07
	Masaoka	ISUM	30.59	26.59
ISSC 2012 Report [27]	Paik	NLFEA	27.34	22.50
		ISUM	25.59	21.97
		Modified P-M	25.67	22.39
		Smith	28.42	22.13
	Wang	NLFEA	31.00	25.00
		Smith	29.85	25.01
	UoG	RINA Rules	28.20	21.70
Present		Method 1	29.20	21.90
		Method 2	29.20	21.90
		Method 3	29.20	21.90

, and the results for CASE 2 and CASE 3 are shown in

Table 6. Residual strength of CASE 2.

Methodology	Residual Strength–Hogging (MN·m)	Residual Strength–Sagging(MN·m)
Method 1	25.53	17.80
Method 2	23.10	17.89
Method 3	23.08	17.89

and

Table 7. Residual strength of CASE 3.

Methodology	Residual Strength–Hogging (MN·m)	Residual Strength–Sagging(MN·m)
Method 1	25.82	24.64
Method 2	25.65	23.62
Method 3	25.64	23.62

.

It can be seen that, for CASE 1, the results obtained by the three methods for the bending moment and position of the neutral axis are very close. However, there are some differences between the results for CASE 2 and CASE 3. Particularly, the results obtained with Method 1 bear some differences with those obtained with the other two methods, and among them Method 3 is most accurate. For this reason, Method 3 is used to investigate the ultimate limit state function in this study.

4.3. Calculation of the Ultimate Limit State Function

4.3.1. Calculation of Sample Points

Since it is not possible to accurately predict the conditions of a ship at sea, a fit to the ultimate limit state curve is to be made for the ship under different loading conditions. For which purpose, a number of samples of ultimate vertical bending moments and horizontal bending moments need to be obtained. In this study, the ultimate/residual strengths of the ship under different combined loads are obtained for various heeling angles, and then decomposed into the vertical and horizontal bending moments.

The double hull VLCC mentioned before is used for the investigation. The hull is inclined from 0 to 170° with an interval of 10 degrees in order to guarantee the accuracy without incurring too much computation load. The ultimate/residual hogging and sagging strengths of the intact hull and also the damaged hull due to collision or stranding are calculated for each angle, and the sample points of the ultimate vertical and horizontal bending moment under combined load are obtained. For each case of intact/damaged hull, 36 sample points are obtained, as shown in Figure 17.

It can be seen in Figure 17 that, in the case of an intact hull, the ultimate bending moment envelope curve is smooth and approximately symmetric; however, in the cases of a damaged hull, there are abrupt changes in the distribution of the samples, and the curves are not symmetric, particularly in the second quadrant where the horizontal bending moment acting on the inclined ship is even larger than the horizontal bending moment under pure horizontal load.

4.3.2. Fitting Results

The fitted curves which obtained the three methods are shown in Figure 18, Figure 19 and Figure 20, where the results are nondimensionalised by dividing them by the pure bending moment. In these figures, Sample point represents the nondimensional results, Fitting 1 the explicit method, Fitting 2 the polynomial fitting method, and Fitting 3 the improved polynomial fitting method. Since the values of the sample points must be positive in the explicit methods, the absolute values are used here. Because the samples points are distributed in all the four quadrants, and they cannot be fitted with these fitting methods using one function, the fitting is made to the sample points located in the same quadrant with the two on the adjacent axes (10 sample points in total). Thus, four fitting functions are obtained for each case.

The results obtained with the three methods for the fitting for the sample points are shown in Figure 21, Figure 22 and Figure 23. In order to compare the accuracy of the fitting in the ranges outside the sample points, the sample points are removed and the results are shown in Figure 24, Figure 25 and Figure 26.

In order to compare the accuracy of the fitting methods at the sample points, the actual vertical bending moments and the fits to them are plotted in Figure 21, Figure 22 and Figure 23, where the abscissa is the fit to the vertical bending moment, and the coordinate is the actual bending moment. If the point is located on the $y=x$ line, it means that the fit is equal to the actual bending moment; and if the larger is the distance between the point and $y=x$, the larger is the difference between the fit and the actual bending moment.

To compare the accuracy of the fitting in the ranges outside the sample points, the fitting is made with the sample point of concern removed; i.e., the sample point is removed, and the fitting is made with the remaining sample points. The vertical bending moment is obtained by substituting the horizontal bending moment of this sample point into the corresponding fitting function, and then compared with the actual vertical bending moment. The same is done for all the sample points, and the comparisons are shown in Figure 24, Figure 25 and Figure 26, where the abscissa is the vertical bending moment calculated from the fitting function, and the coordinate is the actual moment.

In order to evaluate the goodness of fit obtained by each method, the coefficient of determination $R^2$ is calculated using

$$R^2=1-\frac{{\displaystyle \sum {(y-\widehat{y})}^2}}{{\displaystyle \sum {(y-\overline{y})}^2}}$$

(18)

where $y$ is the fit, $\widehat{y}$ the actual value, and $\overline{y}$ the mean of $\widehat{y}$. The maximum of $R^2$ is 1, and the closer $R^2$ is to 1, the better is the fit. The results obtained for the three cases are shown in

Table 8. Goodness of fit for the case of the intact hull.

Fitting No.	1st Quadrant		2nd Quadrant		3rd Quadrant		4th Quadrant
	${\mathit{R}}^2$	Satisfactory?	${\mathit{R}}^2$	Satisfactory?	${\mathit{R}}^2$	Satisfactory?	${\mathit{R}}^2$	Satisfactory?
Fitting 1	1.000	YES	0.999	YES	0.999	YES	1.000	YES
Fitting 2	0.999	YES	0.999	YES	0.999	YES	1.000	YES
Fitting 3	1.000	YES	1.000	YES	1.000	YES	1.000	YES

,

Table 9. Goodness of fit for the case of the damaged hull due to collision.

Fitting No.	1st Quadrant		2nd Quadrant		3rd Quadrant		4th Quadrant
	${\mathit{R}}^2$	Satisfactory?	${\mathit{R}}^2$	Satisfactory?	${\mathit{R}}^2$	Satisfactory?	${\mathit{R}}^2$	Satisfactory?
Fitting 1	1.000	YES	0.981	NO	1.000	YES	0.986	NO
Fitting 2	0.999	YES	0.983	NO	1.000	YES	0.974	NO
Fitting 3	1.000	YES	1.000	YES	1.000	YES	1.000	YES

and

Table 10. Goodness of fit for the case of the damaged hull due to stranding.

Fitting No.	1st Quadrant		2nd Quadrant		3rd Quadrant		4th Quadrant
	${\mathit{R}}^2$	Satisfactory?	${\mathit{R}}^2$	Satisfactory?	${\mathit{R}}^2$	Satisfactory?	${\mathit{R}}^2$	Satisfactory?
Fitting 1	0.998	NO	1.000	YES	0.997	NO	0.999	YES
Fitting 2	0.997	NO	1.000	YES	0.997	NO	0.999	YES
Fitting 3	1.000	YES	1.000	YES	1.000	YES	1.000	YES

. In this study, a fitting method is considered as satisfactory if $R^2\ge 0.999$. The results obtained for the case are given in Table 8, Table 9 and Table 10.

4.4. Discussions

It can be seen that the bending moment is symmetric for the intact ship, but asymmetric for damaged cases. The bending moment of the damage hull is smaller than that of the intact ship; however, some sample points of the damaged case are located outside the bending moment envelope of the intact ship. Fitting for the nondimensionalised bending moment shows that the three fitting methods produce similar results, except that some differences appear when the load is (almost) a pure bending moment. The fitting results are basically those symmetrical about the axes, which means that a symmetry exists for the bending moment in the horizontal plane as well as in the vertical plane (i.e., hogging and sagging bending moments are symmetric). In the cases of the damaged hull, there are some differences between the three methods, particularly, in the results obtained by the classic method and the polynomial fitting method for the second quadrant.

For the case of an intact hull, the bending moment sample points are symmetrical as the cross-section is symmetrical about the centre plane. It can be seen that the bending moments are larger than that of the damaged hull. The sample points of the damaged hull are asymmetric due to the asymmetric cross-section. Especially in the case of a damaged hull due to collision, where the damage is close to the side wall, a clear asymmetric pattern is displayed. Because of the redistribution of stress in the damaged cross-section, some of the bending moments are relatively large and located outside the fitted curve, as shown in Figure 19. Nondimensionalisation of the ultimate bending moments shows that the fitted curve of the samples points of the intact hull is symmetric about the abscissa and also the ordinate. The distribution of the sample points of the stranding case is somewhat asymmetric, but still similar to the case of the intact ship, and this is because the damage is close to the middle of the bottom, and its influence on the vertical bending is limited. However, for the horizontal bending moment, the influence is much more distinct. For the case of collision, the damage is located at the sidewall and part of the deck is missing, and absence of these strong structure members will result in large reduction of the load carrying capacity, and in turn the pattern of distribution is different.

For sample points with small variations, e.g., in the case of the intact hull, all the fitting methods are accurate. However, when there are large variations in the sample points, e.g., in the case of a damaged hull due to collision, the accuracy of Fitting 1 and Fitting 2 is poor, while that of Fitting 3 is still satisfactory. As can be seen in Figure 21, the accuracy of fitting with the explicit method and the polynomial fitting method is poor in the second and fourth quadrant.

Comparison of Figure 21, Figure 22 and Figure 23 with Figure 24, Figure 25 and Figure 26 shows that the fitted points obtained with the corresponding sample points removed are located at larger distances from $x=y$, and Fitting 3 is always more accurate than Fitting 1 and Fitting 2.

Comparison of the results obtained with and without removing the sample points shows that increase in the number of sample points improves the accuracy, but the improvement is limited in many cases; in particularly, in the case of a damaged hull due to collision.

The results obtained for the goodness of fit also shows that the accuracy of Fitting 1 and Fitting 2 is not always accurate for the case of damaged hull due to collision or stranding, while Fitting 3 is the most accurate in all of the cases.

5. Conclusions

In this study, a PSO-based Smith Method is applied to calculate the ultimate strength of the intact/damaged hull under the combined load of the vertical and horizontal bending moment. A weighted piecewise polynomial fitting method is devised to obtain the ultimate limit state function of ships under the combined load of vertical and horizontal bending moment. Based on the comparison of the results with literature, the following conclusions are drawn.

(1)

Compared to the implementations of the Smith Method, where only the translation of the neutral axis is considered, or the linear search method is applied to trace the position of the neutral axis, the PSO-based Smith Method is more accurate in obtaining the instantaneous position of the neutral axis, and helps produce more accurate results for the ultimate bending moment. Application of the PSO-based method also helps with more accurate fitting of the ultimate limit state functions.

(2)

The weighted piecewise fitting method proposed in this study is more accurate for the fitting of the ultimate limit state function than the explicit method and the original polynomial fitting method, especially in the case of a damaged hull due to collision.

(3)

Comparison of the fitting with and without the sample points removed shows that increasing the number of sample points may help with the accuracy of the explicit fitting and the original polynomial fitting method, but the improvement is limited, and the accuracy is still unsatisfactory in some cases.

(4)

Results obtained for the cases of the intact/damaged hull show that the explicit fitting method and the original polynomial fitting method produce fairly accurate results for the case of an intact hull. Thus, they can be used as the first choices for obtaining the ultimate limit state function of intact hulls for their simplicity, while the weighted piecewise fitting method can be used for the cases of damaged hulls.