Multi-Crack Interaction and Influence on the Spherical Pressure Hull for a Deep-

Rujun Li · Yongmei Zhu · Wenjing Fang · Baoji Yin

Abstract This study investigates the interaction and influence of surface cracks on the spherical pressure hull of a deep-sea manned submersible. The finite element model of the spherical hull is established, and a semi-elliptical surface crack is inserted in the welding toe of the spherical hull as the main crack. Considering the combined effect of external uniform pressure and welding residual stress at the weld toe, the stress intensity factor (SIF) is obtained based on the M-integral method. Inserting disturbing cracks at different positions on the spherical hull surface, the interaction and influence between multi-cracks are revealed by numerical calculation. The results show that the existence of the disturbing crack has a great influence on the stress intensity factor of the main crack, and the influence is different with the different location of disturbing crack. The study of the interaction of multiple cracks under different interference factors and the influence of disturbing cracks on the main crack can provide some reference for future engineering applications.

Keywords Spherical pressure hull · Surface cracks · Multi-crack interaction · Stress intensity factor · Welding residual stress · Manned submersible · Finite element model

1 Introduction

A spherical pressure hull with a medium thickness is adopted for the deep-sea manned submersible. It is welded by a hemisphere or several spherical parts (Figure 1). Surface cracks at the weld toes are the most common defects in welded structures. When the manned submersible is under the combined action of the seawater external pressure load and the welding residual stress at the welding seam, fatigue cracks are prone to initiate and propagate at the welding toes,resulting in structural failure. In practical structures, multicracks often exist at the same time, and their size, direction, and distribution are arbitrary. Moreover, there will be a shielding or intensification effect on the crack growth trend because of the cracks’ different locations and directions (Fan et al. 2015; Lam and Phua 1991). It is necessary to study the interaction of multi-cracks and its influence law. Currently,the fracture problem with a multi-cracks structure is a key issue in structural integrity research. The stress intensity factor (SIF), as an important parameter to describe the stress field intensity of crack tip, is the primary problem to be solved in the research.

Figure 1 Spherical pressure hull of manned submersible, adopted by Li WY (Li et al. 2016)

At present, the calculation methods of SIF include analytical methods, numerical methods, and experimental methods. Numerical methods are the most commonly used method for solving the SIF of complex three-dimensional fractures. With the improvement of computer technology,the finite element method (FEM) has the advantages of not being restricted by the geometric form of the cracked body and load, wide application range, and high accuracy. Commonly used finite element methods include the M-integral method (Freund 1978), displacement correlation method(Gerstle and Abdalla 1990), and the virtual crack closure technique (VCCT) (Rybicki and Kanninen 1977). Several scholars have focused on the SIF of two-dimensional and three-dimensional cracks under a complex stress field (Li and Wang 2000; Ayhan and Yücel 2011; Stepanova and Roslyakov 2016; Pasca et al. 2013; Xiao and Yan 2007). For example, Lv et al. (2018) evaluated the SIFs for parallel surface cracks on the plate as a function of crack front position,relative position, aspect ratio, and relative size by numerical analysis based on the M-integral method. Perl et al. (Perl et al. 2016; Perl and Steiner 2015, 2018) calculated the 3-D SIFs of internal surface multiple cracks of spherical pressure vessels by FEM. Jie et al. (2017) used the M-integral method to analyze the surface crack at a cruciform welded joint, verifying that this method is correct and efficient to calculate the SIFs. Ai et al. (2018) calculated the SIF distribution at the front edge of a semi-ellipse surface crack on a flat specimen based on the interaction integral method. The relative error between the result and the theoretical solution was within 4%. Wang and Gao (2019) calculated the SIFs for sliding fatigue-wear with the FEM and carried out the crack propagation and wear analysis well. Xiong et al.(2020) took the finite element numerical solution of the SIF at the midpoint of the crack front edge as the input value and used the ABAQUS and FRANC3D joint simulation to calculate the crack growth step size and the fatigue life of the flat specimen. Most studies on 3-D cracks are concentrated on simple components with a single crack, while less attention has been paid to the interaction and influence of multicracks on the spherical pressure hull surface. The spherical pressure hull has several large openings and is welded by the hemisphere equatorial. Moreover, the welding residual tensile stress makes it more susceptible to multiple radial cracking near the welds.

In the present study, the finite element models of the spherical hull are established, and the semi-elliptical surface cracks are inserted at the weld toes of the spherical hull. The SIFs of the surface cracks are studied numerically under external uniform pressure and welding residual stress.By inserting the disturbing cracks at different locations, the influence on the main cracks is analyzed and the interaction of multi-surface cracks under different SIFs is revealed in detail. This paper aims to provide a rational approach to predict the crack growth behavior and fatigue life of spherical pressure hulls.

2 The Analytical and Numerical Results of SIF

Newman and Raju (1981) proposed an empirical formula for calculating the SIF of a surface crack on a flat plate under a tensile force (Eq. (1)). The accuracy of the SIF calculated by this formula is high, but it is only applicable to flat plate structures under a tensile bending load.

whereKIrefers to the mode I SIF,σtrefers to the tensile(compression) stress,σbrefers to the bending stress,arefers to the crack depth,crefers to the half-length of the crack,trefers to the thickness of the plate, andwrefers to the width of the plate.

This study uses a flat plate model (Figure 2) with a central semi-elliptical surface crack to validate the accuracy of the numerical approach. The length, width, and thickness of the plate are 800, 400, and 20 mm, respectively. The surface crack is located at the center of the plate. The crack depthais 0.3, 0.5, 1, 2, and 4 mm, and the ratio of the crack depth and half-lengtha/cis 0.6. The material selected is the 20MnMoNb alloy steel. The parameters are set as follows:yield stress of 490 MPa, Young’s modulus of 210 GPa, and Poisson’s ratio of 0.3. During the calculation, a fully constrained boundary condition is applied at one end of the plate model, and a tensile force of 100 MPa is applied at the other end (Zhu et al. 2020). Surface cracks are introduced and the mesh will then be re-divided as shown in Figure 2b. The mesh of the plate’s finite element model is divided freely and uniformly. The total number of elements is approximately 800, and the number of nodes is approximately 1700. In addition, the 20-node singular elements are divided around the crack front with an element length of 0.01 mm.

Figure 2 Flat plate model with a central surface crack (Zhu et al.2020)

The SIFs of the crack front are solved by the M-integral method (Budiansky and Rice 1973; Knowles and Sternberg 1972; Rice 1964; Wawrzynek et al. 2005) in FRANC3D(Fracture Analysis Consultants Inc. 2016). The finite element numerical solutions are compared with the theoretical solution of the Newman-Raju formula. Tables 1 and 2 show the results of the SIF at the crack surface pointKIcand the deepest pointKIa, respectively.

Table 1 Analytical and numerical solutions of the SIF of the surface point KIc

Table 2 Analytical and numerical solutions of the SIF of the deepest point KIa

The error between the M-integral and the Newman-Raju solution is within 3%. Therefore, the results obtained by the M-integral method for calculating the SIF are accurate and reliable, which is consistent with Ai’s finding (Ai et al. 2018). Thus, the M-integral method is used to calculate the SIF of the surface crack of the spherical pressure hull, and the FRANC3D can be used to do the following analysis.

3 Solutions of SIF

3.1 Finite Element Model

This study introduces a semi-elliptical surface crack near the equator of the spherical pressure hull as a model for finite element analysis. Figure 3 shows the modeling process. First, a finite element model of the spherical pressure hull is established through the ABAQUS, which is then divided into two parts: the sub-model and the remaining model. The sub-model is used to introduce the initial crack. The combination of the sub-model and the remaining model is the global model, which reflects the boundary conditions and load distribution. The established crack node data file is imported into FRANC3D, and each node of the crack front edge is selected as a set, which is convenient to divide the singular element around the crack tip.After the cracks are imported, the sub-model elements are re-divided, preparing for the next step of the finite element calculation process.

Figure 3 Numerical model for the spherical pressure hull with cracks

3.1.1 Spherical Pressure Hull

This study concerns the surface cracks at the equatorial weld toe of a 7000-m manned spherical pressure hull. The spherical hull is made of high-strength titanium alloy. Table 3 lists the geometric parameters and material properties of the model, and Figure 4 shows the simplified numerical model used in the study. For a medium-thickness spherical hull,numerical solutions tend to converge when the ratio of the volume element size to the radius is 0.06-0.07 (Zhang et al.2015). Considering the calculation accuracy and efficiency,the elements are divided into the shape of a tennis ball and the element type is a three-dimensional 8-node reduced integral volume element. The average size of the circumferential element is 0.06R, which is 67.5 mm, and the thickness direction is four layers. The volume element numerical scheme of the model is established with 15,392 elements and 19,250 nodes (Figure 4b).

Table 3 Geometric and material properties of the model

Figure 4 Numerical model of the spherical hull

In addition to the uniform seawater pressure load,the spherical pressure shell is not subjected to any other restrictions when working underwater. The rigid displacement of the structure needs to be eliminated when the finite element model is calculated. Therefore, it is necessary to add corresponding boundary constraints to the model. In this study, six displacement components are constrained by three-point constraints. As shown in Figure 4b, nodes 1 and 2 are taken on theX-axis to limit the displacement along theY- andZ-axis (Uy=Uz=0). Node 3 is taken at the position 90° apart from the two points in the same longitude to limit the displacement along theX- andY-axis(Ux=Uy=0) (Yu et al. 2005).

Figure 5 Finite element model of the spherical hull with a single crack

3.1.2 Insert Initial Surface Crack

This study describes the surface crack by a semi-elliptic crack and the finite element model of a spherical hull with a single crack is established as shown in Figure 5. The size of the semi-elliptical surface cracks are as follows: depth and half-length ratioa/c=1 and crack deptha=4 mm. In this study, the cracks are inserted into the hemispherical welds of the spherical hull.

Because of ther-1/2singularity of the stress field at the crack tip, the crack tip elements need to be divided into singularity elements. Since the SIF of the crack tip that is calculated by the 20-node singular element model is independent of the meshing near the crack tip and the numerical analysis results are stable and reliable (Chen et al. 2010), the 20-node singular elements were divided around the crack front as shown in Figure 6a. In addition, the crack tip element angleα=45° and the number of layers around the crack tip elementntip=3 (see Figure 6a) reduce the discretization errors by symmetrical distribution. Figure 6b shows the crack surface mesh whereais the crack depth andcis the crack half length. By encrypting the crack front mesh to ensure the accuracy of SIFs, the remaining part is coarser to reduce the number of elements, thereby reducing the computational cost.

Figure 6 Crack surface mesh

3.2 Load Analysis

3.2.1 External Pressure

The working depth of the deep-sea manned submersible studied in this study is 7000 m. The outer surface of the spherical hull is subjected to a uniformly distributed pressure. According to Eq. (2):

wherekis the safety factor,ρwis the density of seawater,gis the acceleration of gravity, andHis the working depth of the sea. Since the pressure hull runs smoothly in operation,kis set as 1.35,ρwis 1.025 g/cm3,gis 9.8 m/s2, andHis 7000 m. The calculated load is 105.47 MPa.

Therefore, a load of 105.47 MPa is applied to the finite element model as shown in Figure 4a.

3.2.2 Welding Residual Stress

As a welding structure, the spherical pressure hull will inevitably have welding residual stress at the weld seam. Surface cracks easily occur at the weld toe of the opening, most of which are parallel to the weld direction and are subjected to welding residual stress perpendicular to the crack surface(Al-Mukhtar 2013). For the surface cracks parallel to the weld fusion line, the maximum welding residual stress on the surface of the weld toe is 0.3σs(Huang et al. 2011). The distribution of the welding residual stress on the weld is inhomogeneous. Based on the research results of Ikushima and Shibahara (2015), it is assumed that the welding residual stress perpendicular to the circumferential weld direction is linearly distributed along the plate thickness direction.Moreover, the distribution of the welding residual stress at the weld toe of the spherical hull along the thickness direction can be simplified to a linear distribution. Therefore, if the thickness of the plate ist, the distribution of residual stress along the thickness of the plate can be expressed as the following linear distribution (x=0 is the outer surface of the weld toe):

whereσRis the maximum welding residual stress, which is 0.3σs.

According to Eq. (3), the weld residual stress in the outer surface of the weld toe is tensile residual stress, and the maximum value is 0.3σs. It gradually reduces along the thickness direction (the value at 1/4 thickness is 0), and the compressive residual stress gradually increases, reaching a maximum of 0.3σsat 1/2 thickness. The weld residual stress is symmetrically distributed in the mid-surface of the spherical hull (see Figure 7). Figure 7a shows the residual stress distribution, which is applied to the crack surface. As shown in Figure 7b, a local coordinate system is established at the center of the outer surface of the crack. The crack surface is located in theX-Yplane, and thex-axis points to the deepest point of the crack. Based on the local coordinate system, the residual stress that is linearly distributed along thex-axis is established by FRANC3D. After the FEM analysis is finished, three modes of SIFs (KI, KII, KIII) along the crack front are computed by FRANC3D based on the M-integral.

Figure 7 Distribution of the residual stress on the crack surface

3.3 Results and Discussion

Figure 8 presents the distribution of the SIF (KI,KII, andKIII)values along the crack front, whereθis the angle between the crack front and the crack outer surface (see Figure 7a).It can be seen from Figure 8 that theKIvalue of the crack front is much larger than the values ofKIIandKIII.KIwas 662.54 MPa·mm1/2at point A (θ=0°), 665.87 MPa·mm1/2at point B (θ=180°), and 550.99 MPa·mm1/2at the deepest point (θ=90°), while the maximum values ofKIIandKIIIare only 4.22 MPa·mm1/2and 1.71 MPa·mm1/2, respectively.Therefore, the surface cracks of the spherical hull are mainly characterized by mode I fractures. Compared withKI, the variation ofKIIandKIIIvalues are very small withθ, i.e., the influence of mode II and III SIF on the crack growth is very small. Mode I fracture is the most dangerous; therefore, the SIFKIvalues along the surface cracks are discussed below.Figure 9 shows the variation ofKI_SUM,which is the maximumKIunder the combined action of external pressure and welding residual stress andKI_RS, which is the maximumKIunder the action of welding residual stress. It is evident in Figure 9 that the curves of the two conditions are very close and the maximum difference is only 2.2%. Thus, for mode I SIF, the effect of external pressure on the SIF of the crack front is relatively small, and the main influencing factor is the welding residual stress at the welding toe.

Figure 8 Distribution of KI, KII, and KIII along the crack front

Figure 9 Variation of the maximum KI with crack depth

4 Multi‑Crack Interaction Analysis

4.1 Finite Element Model

This work focuses on the interaction between two surface cracks. For convenience, the disturbing crack inserted is the same as the main crack. Figure 10 shows the finite element model of a spherical hull with multiple surface cracks.

In this study, some examples of multi-cracks are analyzed. The variation of mode I SIFs along the main crack fronts and the interaction effect on the crack tipsA,B, and the deepest point (see Figure 10b) of the main crack are discussed after inserting the disturbing crack.

4.2 Disturbing Parameter Definition

The relative positions of the disturbing crack and the main crack on the surface of the spherical hull can be expressed by the position parametersαandβwhereαis the angle in theX-Yplane formed by the line connecting the center of the main crack to the center of the sphere and the other line connecting the center of the disturbing crack to the center of the sphere, andβis the angle in theY-Zplane formed by the line connecting the center of the main crack to the center of the sphere and the other line connecting the center of the disturbing crack to the center of the sphere (Figure 11).

Figure 10 Finite element model of a spherical hull with multi-surface cracks

Figure 11 Position parameter definition

Figure 12 Variation of KI along the main crack fronts with different α

The variation of SIF at the right end pointAof the main crack was observed by adding the disturbing cracks at different positions on the surface of the spherical hull. To analyze the effect of the disturbing crack on the main crack, the SIF reductionλis defined as:

whereK0is the SIF of the main crack without the disturbing crack andKis the SIF of the main crack with the disturbing crack.

4.3 Results and Discussion

Example 1The disturbing crack is inserted whenβ=0°, i.e.,the disturbing crack and the main crack are on the equatorial plane. Since the diameter of the spherical hull is 2250 mm,αis 0.45° when the distance between the two cracks is 9 mm.Thus, the influence of the change ofαonKIalong the main crack fronts is analyzed whenαis 0.45°, 0.5°, 0.6°, 0.8°,and 1°. Figure 12 presents the variation ofKIalong the main crack fronts with differentα. It can be seen from Figure 12 that theKIvalues along the main crack fronts near pointAincrease in varying degrees with the decrease ofα, while theKIvalues near point B hardly change. The increase ofKIis the largest at pointA. Table 4 lists theKIandλvalues of pointA. Whenαis 0.45°,KIis 776.06 MPa·mm1/2, which is 15.6% larger than theKIvalue of the single crack. It can be seen that whenαis greater than or equal to 0.8°, the SIF reductionλis less than 0.007, which means that the disturbing crack has little influence on the main crack. Figure 13 presents the influence curvesλ-αon pointsAandBand the deepest point of the main crack. It can be seen thatλincreases with the decrease ofα, and the increase of range becomes larger at pointA. The SIF intensification of pointBand the deepest point by the disturbing crack is very weak.The influence on pointsA,B, and the deepest point with the variation ofαis fully evaluated by theλ-αcurves.

Example 2The disturbing crack is inserted whenα=0° and is located directly above the main crack. The influence of the disturbing crack onKIalong the main crack fronts is analyzed whenβis 0.02°, 0.05°, 0.1°, 0.3°, 0.5°, and 1°.Figure 14 presents the variation ofKIalong the main crack fronts with differentβvalues. It is evident in Figure 14 that with a decrease inβ, theKIvalues become smaller along the main crack fronts and the shielding effect becomes greater on the two cracks. This means that when the disturbing crack is in this position, crack initiation will be shielded. Moreover, the closer the distance is, the greater is the shielding effect. Table 5 lists theKIand λ values of point A.KIis 403.12 MPa·mm1/2whenβis 0.02°, which is 39.2% less than theKIvalue of the single crack. Figure 15 shows the influence curvesλ-βon pointsA, B, and the deepest point of the main crack. As shown in Figure 15, the curvesλA,λB, andλdeephave the same trend. Whenβis greater than 1,λgradually tends to 0, i.e., the disturbing crack has little influence on the main crack. Therefore, the influence of the disturbing crack at this position can be ignored in the engineering process, and such cracks can even be prefabricated to reduce the SIFs.

Example 3To further analyze the inhibitory action of the disturbing crack on the main crack, the influence of the disturbing crack onKIalong the main crack fronts was analyzed whenβ=0.1° andαwas 0, 0.05°, 0.1°, 0.2°, 0.3°, 0.4°, 0.5°,and 0.8°. Figure 16 presents the variation ofKIalong themain crack fronts with differentαvalues and Table 6 lists theKIand λ of pointA. It can be seen that whenαis less than 0.4°, namely, the distance between pointsBandCis less than 8 mm (i.e., 2c, wherecis the half-length of the main crack), the SIFs of the main crack are less than that of the single crack. This indicates that the growth of the main crack is shielded by the disturbing crack. Whenαis equal to 0.2° (the distance betweenBandCis approximately equal toc), the shielding effect on pointAis the greatest.KIis 203.06 MPa·mm1/2, which is 69.3% smaller than theKIvalue of the single crack. Whenαis larger than 0.4°,i.e., the distance between points B and C is larger than 2c,the SIFs of the main crack are intensified by the disturbing crack. With the continuous increase ofα, the intensification of the disturbing crack on the main crack becomes weaker because of the long distance between the two cracks. From theλ-αcurves in Figure 17, it can be seen that the influence trend of the disturbing crack on the main crack is basically the same at point B and the deepest point, and the influence range is much smaller than that of pointA. According to the conclusion in Example 1, with the increase ofα, the values ofλA,λB, andλdeepwill all tend to 0.

Figure 13 λ-α curves of points A, B, and the deepest point affected by the disturbing crack

Figure 14 Variation of KI along the main crack fronts with different β

Table 4 KI and λ of point A

Example 4In Examples 1 to 3, the main crack and the disturbing crack were parallel. In example 4, the influence of the disturbing crack on the main crack is investigated when the two cracks are non-parallel with a certain angle.A parameterφ(Figure 18) is defined as the angle between two crack surfaces. Here,αis 0.5° andβis 0°. The influence of the disturbing crack on the main crack is analyzed with differentφ. Figure 18 presents the variation ofKIalong the main crack fronts with differentφand Table 7 lists theKIandλof pointA. As shown in Figure 18, theKIvalue of pointA(704.18 MPa·mm1/2) is the largest whenφis 0°,which is 6.4% larger than thatof the single crack. Therefore,the disturbing crack has the strongest enhancement effect on the main crack when the two cracks are coplanar. Figure 19 shows the influence curvesλ-φon pointsA, B, and the deepest point of the main crack. As can be seen,λBandλdeepvalues are less than 0.01 with differentφ, indicating that they are insensitive to the variation ofφ. The curvesλAdecreases with the increase ofφ, but it still shows a weak intensification effect in the middle range ofφbetween 0° and 90°.

Moreover, whenαis 0° andβis 0.3°, the influence of the disturbing crack on the main crack is analyzed with differentφ. Figure 20 presents the variation ofKIalong the main crack fronts with differentφ, and Table 8 lists TheKIand λ values of pointA. As shown in Figure 20, theKIof the main crack is less than that of the single crack whenφis 0° and 30°. This suggests that the disturbing crack shields the main crack initiation. The effect is the strongest whenφis 0°. The disturbing crack enhances at pointAand shields at pointBwhenφis 60°. TheKIvalue of pointAincreases by 6% whenφis 60°. The disturbing crack shields the two crack tips (AandB) of the main crack and intensifies the deepest pointwhenφis 90°. TheKIvalue of the deepest point increases by 3.5% whenφis 90°. Figure 21 shows the influence curvesλ-φon pointsA, B, and the deepest point of the main crack.As shown in Figure 21, theλdeepvalue increases with the increase ofφ, i.e., from-0.22 (KI=429.76 MPa·mm1/2,φ=0°) to 0.036 (KI=570.56 MPa·mm1/2,φ=90°).

Figure 21 λ-φ curves of points A, B, and the deepest point affected by the disturbing crack

Table 8 KI values and λ of point A

Figure 19 λ-φ curves of points A, B, and the deepest point affected by the disturbing crack

Figure 20 Variation of KI along the main crack fronts with different φ (α=0°, β=0.3°)

Table 7 KI values and λ of point A

Figure 17 λ-α curves of points A, B and the deepest point affected by the disturbing crack (β=0.1°)

Figure 18 Variation of KI along the main crack fronts with different φ (α=0.5°, β=0°)

Table 6 KI and λ of point A

Figure 15 λ-β curves of points A, B, and the deepest point affected by the disturbing crack

Figure 16 Variation of KI along the main crack fronts with different α (β=0.1°)

Table 5 KI and λ values of point A

5 Conclusions

(1) The surface cracks of the spherical hull are mainly characterized by mode I cracks under the combined action of seawater pressure and welding residual stress.Compared with mode I SIF, modes II and III SIFs have little effect on the cracks. Moreover, the existence of the disturbing crack has a great influence on the SIFs of the main crack. The different locations of the disturbing crack have a great different influence, which is consistent with practical engineering.

(2) Whenβis 0° andφis approximately 0°, which means that the disturbing crack is coplanar with the main crack,the closer is the distance between the multi-crack tips,the higher is the fracture enhancement. Whenαis 0.5°and the two crack faces have an angleφ, the effect of the enhancement changes with the change ofφ. Whenφis 0°, the disturbing crack has the strongest enhancement effect on the main crack. Thus, the SIF intensification is largest when the two cracks are coplanar, which will greatly accelerate the propagation speed between the cracks before they merge into a larger crack.

(3) Whenαis 0° andφis approximately 0°, the disturbing crack shields the crack tips and the deepest points of the main crack. The closer is the distance between the multicracks, the greater is the shielding effect. Whenβis 0.1°andαchanges, the distance between the ipsilateral tips of two cracks is less than 2c(the length of the main crack)and the shielding effect occurs. The shielding effect on point A is the greatest when the distance between the ipsilateral tips of the two cracks is approximatelyc, and itsKIdecreases by 69.3% than that of the single crack.Therefore, the influence of the disturbing crack at these positions can be ignored in practical engineering, and the SIFs can be reduced by prefabricating such cracks. When the distance between the ipsilateral tips of two cracks is larger than 2c, the enhancement effect occurs. Whenφis less than 30°, the shielding effect occurs. Otherwise, the enhancement effect occurs.

(4) The influence of the disturbing crack on the main crack can guide the crack in engineering practice. The interaction of double cracks is discussed in this study,which awaits further investigation in view of three or more cracks.

FundingThis work was supported by the State Key Program of National Natural Science of China (Grant No. 51439004), the Natural Science Foundation of the Jiangsu Higher Education Institutions of China (Grant No. 19KJA530002), and the Open Project Foundation of State Key Laboratory of Ocean Engineering (No. 2006).