Investigation of the Dynamic Buckling of Spherical Shell Structures Due to Subsea Collisions

: This paper is the ﬁrst to present the dynamic buckling behavior of spherical shell structures colliding with an obstacle block under the sea. The effect of deep water has been considered as a uniform external pressure by simplifying the effect of ﬂuid–structure interaction. The calibrated numerical simulations were carried out via the explicit ﬁnite element package LS-DYNA using different parameters, including thickness, elastic modulus, external pressure, added mass, and velocity. The closed-form analytical formula of the static buckling criteria, including point load and external pressure, has been ﬁrstly established and veriﬁed. In addition, unprecedented parametric analyses of collision show that the dynamic buckling force (peak force), mean force, and dynamic force redistribution (skewness) during collisions are proportional to the velocity, thickness, elastic modulus, and added mass of the spherical shell structure. These linear relationships are independent of other parameters. Furthermore, it can be found that the max force during the collision is about 2.1 times that of the static buckling force calculated from the analytical formula. These novel insights can help structural engineers and designers determine whether buckling will happen in the application of submarines, subsea exploration, underwater domes, etc.


Introduction
Deep sea submarines have been used widely in many applications such as investigation, the exploration of oceanographic resources, military action, and so on [1][2][3][4][5]. In reality, all subsea structures are subjected to significant hydrostatic pressure [6,7]. Therefore, the pressure hull is the most important component of any deep submersible. It cannot only offer the capability for resisting the external load, it can also provide more workspace for captains or scientists. The shape of the hull is usually designed as a spherical shape [8,9]. In practice, the exploration vehicle usually works at a sea depth of 1.5-8.5 km under the mean sea level with a low speed [10][11][12][13]. During operations, a collision could occur between the submarine and any obstacle such as a rock wall, iceberg, cliff, or another submarine, as shown in Figure 1. The collisions can induce the damage of the submarine nose structure through the failure of components, and the unstable deformation or buckling of the spherical nose. The buckling behavior and buckling criterion of spherical shell structures subjected to collisions have attracted more research attention over recent years [14][15][16].
Many previous studies have been carried out that investigate spherical shells under a uniform pressure [17][18][19][20]. The first and most preeminent researcher who established the buckling theory is Timoshenko [21]. His classical conclusion was derived under the elastic material and small deformation deformation assumption. Afterwards, many studies considering large deformation and more complicated loads have been investigated [22][23][24]. Smith [25] evaluated the buckling of a multisegment underwater hull. In his study, it was revealed that the collapse pressures from numerical analyses agreed well with all five experimental results obtained from the collapse test of a laboratoryscale mild steel structure. Zhang and Wang [26] studied the effect of thickness on the buckling strength of egg-shaped pressure hulls. The results showed that the egg-shaped pressure hulls seem to be applicable to deep sea manned/unmanned submersibles, especially to full ocean depth ones (11,000 m). Khakina [27] and Yu [28] studied the buckling load using the energy method for subsea collisions. However, they studied the static problem without dynamic effect. It was reported that the deep pressure hull could be buckled without exceeding the limited strength under collision [29,30], which can be strongly affected by the shell thickness, collision velocity, and properties of materials [31,32]. Nevertheless, spherical pressure hulls are the most efficient and popular type for the deep-sea pressure hull structures [33][34][35], as the spherical hull has the lowest buoyancy factor (weight-to-buoyancy ratio) [36,37]. Liang and Shieh [38] reported the results of a numerical study on the optimal design of a multi-segment spherical hull, which contained details regarding the sensitivity analysis of the design variables, including the thickness and the inner radius of the rib-ring. Furthermore, Xu bai et al. [39] analyzed the capacity of the ring-stiffened cylindrical shell via the experimental method. The paper exhibited the influence on the stability of the cylindrical shell, but had little analytical discussion based on the derivations of the buckling loads of the cylindrical shell. Mackay et al. [40] studied the stability of submarines with and without artificial failure. Several numerical models to evaluate stability were established, but whether or not the results could be extended to the other applications was not discussed. Bich et al. [41] analyzed the nonlinear static and buckling of spherical shells under external pressure, incorporating the effects of temperature. They took von Karman-Donnell geometrical nonlinearity and the external pressure load without a concentrated load into account.
Based on the critical literature review, it appears that very few studies exist for the spherical shells that consider both external pressure and collision [42][43][44]. Therefore, this present paper is the first to reveal the dynamic buckling behavior for spherical shell structures experiencing the hydrostatic pressure and the concentrated dynamic force induced by subsea collision. This study employs ANSYS18.0 and LS-DYNA to establish the static and dynamic models, respectively. Parametric analysis has been carried out to determine and validate an explicit formulation for static critical buckling load. To obtain the relationship between dynamic buckling force and the static buckling force, the max and mean dynamic force over the duration of the collision are analyzed. The effects of the material and geometrical properties, external pressure, and added mass of the spherical shells are then investigated and discussed. It was reported that the deep pressure hull could be buckled without exceeding the limited strength under collision [29,30], which can be strongly affected by the shell thickness, collision velocity, and properties of materials [31,32]. Nevertheless, spherical pressure hulls are the most efficient and popular type for the deep-sea pressure hull structures [33][34][35], as the spherical hull has the lowest buoyancy factor (weight-to-buoyancy ratio) [36,37]. Liang and Shieh [38] reported the results of a numerical study on the optimal design of a multi-segment spherical hull, which contained details regarding the sensitivity analysis of the design variables, including the thickness and the inner radius of the rib-ring. Furthermore, Xu bai et al. [39] analyzed the capacity of the ring-stiffened cylindrical shell via the experimental method. The paper exhibited the influence on the stability of the cylindrical shell, but had little analytical discussion based on the derivations of the buckling loads of the cylindrical shell. Mackay et al. [40] studied the stability of submarines with and without artificial failure. Several numerical models to evaluate stability were established, but whether or not the results could be extended to the other applications was not discussed. Bich et al. [41] analyzed the nonlinear static and buckling of spherical shells under external pressure, incorporating the effects of temperature. They took von Karman-Donnell geometrical nonlinearity and the external pressure load without a concentrated load into account.
Based on the critical literature review, it appears that very few studies exist for the spherical shells that consider both external pressure and collision [42][43][44]. Therefore, this present paper is the first to reveal the dynamic buckling behavior for spherical shell structures experiencing the hydrostatic pressure and the concentrated dynamic force induced by subsea collision. This study employs ANSYS18.0 and LS-DYNA to establish the static and dynamic models, respectively. Parametric analysis has been carried out to determine and validate an explicit formulation for static critical buckling load. To obtain the relationship between dynamic buckling force and the static buckling force, the max and mean dynamic force over the duration of the collision are analyzed. The effects of the material and geometrical properties, external pressure, and added mass of the spherical shells are then investigated and discussed.

Static Buckling Criteria
From the introduction above, the buckling critical pressure of the spherical shell under a uniform external pressure load based only on the elastic material can be written as Equation (1) [45]: where E, υ are the elastic modulus and the Poisson ratio of the materials, respectively. t and R are the thickness and the radius of the shell, respectively. However, surprisingly, it is important to note that there is no buckling critical force for the spherical shell under concentrated load only.
By defining p* in Equation (2), where p is the concentrated load, it can be found [45] that the slope of reaction force-deflection curve changed rapidly at around p* = 2.0. Thus, there is no analytical buckling critical load for the spherical shell under both external pressure and concentrated load [46][47][48][49]. From the report [45] mentioned above, the critical force p c as a critical point load situation can be defined by p* = 2.2. Clearly, the critical point load varies with different pressure, so that the buckling criteria is the function of both point load p and external uniform pressure q. In this paper, the critical formulation can be expressed as Equation (3), satisfying two limitations: Considering that there is no buckling when only the shell is subjected to the point load, it should be noted that the N should be a very large or infinite value. The function can be determined by experimental or numerical simulation data via regression analysis [50,51].

Verification of FEM Model
The evaluation result of the spherical shell under concentrated load is shown in Figure 2, which can be verified by the FEM model. The model is based on a perfect hemispherical shell with a radius of 5 m and a thickness of 0.1 m. The displacement has been applied at the apex of the shell, with the bottom of the shell fixed in all of the X-Y-Z directions, as shown in Figure 2a. The ideal elastic-plastic model has been adopted as the constitutive model in which the elastic modulus is 210 GPa, the yield stress is 750 MPa, and the Poisson ratio is 0.23, as shown in Figure 2b. The Shell181 element in ANSYS 18.0 is used with the large displacement option on. Figure 2c shows a comparison plot of the reaction force of the apex versus the deflection ratio, where the deflection ratio was defined as ω/t, and p* was defined as shown in Equation (2).

Static Buckling Criteria
From the introduction above, the buckling critical pressure of the spherical shell under a uniform external pressure load based only on the elastic material can be written as Equation (1) [45]: where E ,  are the elastic modulus and the Poisson ratio of the materials, respectively. t and R are the thickness and the radius of the shell, respectively. However, surprisingly, it is important to note that there is no buckling critical force for the spherical shell under concentrated load only.
By defining p* in Equation (2), where p is the concentrated load, it can be found [45] that the slope of reaction force-deflection curve changed rapidly at around p* = 2.0. Thus, there is no analytical buckling critical load for the spherical shell under both external pressure and concentrated load [46][47][48][49]. From the report [45] mentioned above, the critical force c p as a critical point load situation can be defined by p* = 2.2. Clearly, the critical point load varies with different pressure, so that the buckling criteria is the function of both point load p and external uniform pressure q . In this paper, the critical formulation can be expressed as Equation (3), satisfying two limitations: Considering that there is no buckling when only the shell is subjected to the point load, it should be noted that the N should be a very large or infinite value. The function can be determined by experimental or numerical simulation data via regression analysis [50,51].

Verification of FEM Model
The evaluation result of the spherical shell under concentrated load is shown in Figure 2, which can be verified by the FEM model. The model is based on a perfect hemispherical shell with a radius of 5 m and a thickness of 0.1 m. The displacement has been applied at the apex of the shell, with the bottom of the shell fixed in all of the X-Y-Z directions, as shown in Figure 2a. The ideal elastic-plastic model has been adopted as the constitutive model in which the elastic modulus is 210 GPa, the yield stress is 750 MPa, and the Poisson ratio is 0.23, as shown in Figure 2b. The Shell181 element in ANSYS 18.0 is used with the large displacement option on. Figure 2c shows a comparison plot of the reaction force of the apex versus the deflection ratio, where the deflection ratio was defined as / t  , and p* was defined as shown in Equation (2).  Figure 2c shows the comparison of the reaction force at the apex point versus the corresponding displacement between the numerical data and experimental results from Sabir [49]. In order to demonstrate the influence of the discretization, the numerical result with different meshing size is displayed. Since the results of size = 0.3 m and size = 0.1 m are the same, the size 0.3 m has thus been adopted in this paper. It can be concluded that the results agree very well, and the result trend is not linear or bilinear, but rather is irregular. At the beginning, the slope is steep, and after p' at around 2.2, the slope becomes very mild. Thus, the value of c p can be defined by p* = 2.2.

Model
Several numerical simulations have been conducted to validate the precise expression of buckling formulation and the relevant parameters by regression analysis. The models have simulated identical parameters to the initial studies with thicknesses of 0.05 m and 0.10 m, respectively. In order to obtain the critical load data, including pressure and point load, the different uniform pressure and displacement have been applied 1.0 m at the apex [52][53][54]. The bottom of the shell is set to be fixed, as shown in Figure 2a. Figure 3 shows the typical results of reaction force versus the corresponding displacement. When the displacement is at 0.10 m, the reaction force reaches its maximum value. Then, it decreases with the increment of displacement, implying that the shell has been buckled. Note that the point load has to be in the inverse direction (negative value) in order to enable the large deformation. The buckling load can be obtained as p = 9.8 MN with q = 16 MPa. In addition, the von Mises stress at the apex versus the displacement is depicted in the same figure. The maximum stress is about 415 MPa,  Figure 2c shows the comparison of the reaction force at the apex point versus the corresponding displacement between the numerical data and experimental results from Sabir [49]. In order to demonstrate the influence of the discretization, the numerical result with different meshing size is displayed. Since the results of size = 0.3 m and size = 0.1 m are the same, the size 0.3 m has thus been adopted in this paper. It can be concluded that the results agree very well, and the result trend is not linear or bilinear, but rather is irregular. At the beginning, the slope is steep, and after p' at around 2.2, the slope becomes very mild. Thus, the value of p c can be defined by p* = 2.2.

Model
Several numerical simulations have been conducted to validate the precise expression of buckling formulation and the relevant parameters by regression analysis. The models have simulated identical parameters to the initial studies with thicknesses of 0.05 m and 0.10 m, respectively. In order to obtain the critical load data, including pressure and point load, the different uniform pressure and displacement have been applied 1.0 m at the apex [52][53][54]. The bottom of the shell is set to be fixed, as shown in Figure 2a. Figure 3 shows the typical results of reaction force versus the corresponding displacement. When the displacement is at 0.10 m, the reaction force reaches its maximum value. Then, it decreases with the increment of displacement, implying that the shell has been buckled. Note that the point load has to be in the inverse direction (negative value) in order to enable the large deformation. The buckling load can be obtained as p = 9.8 MN with q = 16 MPa. In addition, the von Mises stress at the apex versus the displacement is depicted in the same figure. The maximum stress is about 415 MPa, which is less than the yield stress, considering that the structure may be failed for either exceeding the limited strength or buckling without exceeding the limited strength, so that the material behaves as an elastic model before buckling. which is less than the yield stress, considering that the structure may be failed for either exceeding the limited strength or buckling without exceeding the limited strength, so that the material behaves as an elastic model before buckling.

Results
According to the model described above, the buckling forces with corresponding pressures in accordance with the different parameters can be obtained. All of these simulation results with different uniform pressure and concentrated load at the beginning of buckling are shown in Tables 1  and 2.

Results
According to the model described above, the buckling forces with corresponding pressures in accordance with the different parameters can be obtained. All of these simulation results with different uniform pressure and concentrated load at the beginning of buckling are shown in Tables 1 and 2.

Buckling Criteria Analysis
All of the results of regression analyses are shown in Table 3. The average values of α and β are 0.208 and 0.875. For convenience in application, the parameters of Equation (5) can be admitted as α = 0.2 and β = −0.9.

Discussion on Static Buckling
According to the results of Table 3, Equation (4) can be re-written as Equation (5). It should be noted that when the force p and pressure q satisfy the formula in Equation (5), the static buckling will incur. The pressure p tends to be an infinite value, as the pressure q tends to be zero. In such a case, the buckling force does not exist if the pressure is very small. Similarly, Sabir [49] also showed that when ' q < 0.05, the buckling would not be induced.

Buckling Criteria Analysis
All of the results of regression analyses are shown in Table 3. The average values of α and β are 0.208 and 0.875. For convenience in application, the parameters of Equation (5) can be admitted as α = 0.2 and β = −0.9.

Discussion on Static Buckling
According to the results of Table 3, Equation (4) can be re-written as Equation (5). It should be noted that when the force p and pressure q satisfy the formula in Equation (5), the static buckling will incur. The pressure p tends to be an infinite value, as the pressure q tends to be zero. In such a case, the buckling force does not exist if the pressure is very small. Similarly, Sabir [49] also showed that when q < 0.05, the buckling would not be induced.

Model
The dynamic model has been established by using two segments, including a spherical shell and a cylinder, as depicted as Figure 5a. The nose structure of the submarine can be simplified as a model of a spherical shell with added mass, as illustrated as Figure 5b. The obstacle block in both models has been simulated by solid elements. The shell structure of the submarine nose is modeled using shell elements in LS-DYNA [55]. The material for the obstacle block is set to be the rigid material with E = 35 GPa and a density of 2400 kg/m 3 . For the cylinder, it is usually reinforced by ring bars, a stiffener, or other measures. It is also considered as a rigid body with a thickness of 0.04 m and a density of 7800 kg/m 3 . The spherical shell in both models is assumed to behave as an elastic material with E = 210 GPa and a density of 7800 kg/m 3 . The added mass is set to 35,800 kg, which can be calculated from the total mass of the cylinder model. In addition, the effect of the strain rate was neglected in this study on the basis of the low velocity impacts. In the modeling, the radius of the shell is set to be 2.0 m, with a thickness of 0.03 m and an external pressure of 20 MPa, which implies that the submarine is considered to be in operations at the depth of 2.0 km under the sea.

Model
The dynamic model has been established by using two segments, including a spherical shell and a cylinder, as depicted as Figure 5a. The nose structure of the submarine can be simplified as a model of a spherical shell with added mass, as illustrated as Figure 5b. The obstacle block in both models has been simulated by solid elements. The shell structure of the submarine nose is modeled using shell elements in LS-DYNA [55]. The material for the obstacle block is set to be the rigid material with E = 35 GPa and a density of 2400 kg/m 3 . For the cylinder, it is usually reinforced by ring bars, a stiffener, or other measures. It is also considered as a rigid body with a thickness of 0.04 m and a density of 7800 kg/m 3 . The spherical shell in both models is assumed to behave as an elastic material with E = 210 GPa and a density of 7800 kg/m 3 . The added mass is set to 35,800 kg, which can be calculated from the total mass of the cylinder model. In addition, the effect of the strain rate was neglected in this study on the basis of the low velocity impacts. In the modeling, the radius of the shell is set to be 2.0 m, with a thickness of 0.03 m and an external pressure of 20 MPa, which implies that the submarine is considered to be in operations at the depth of 2.0 km under the sea. The external pressure can be calculated by Equation (6): where  , g , and h are the density of water and gravity acceleration and the height under the sea, respectively. The pressure is then applied to the surface of the shell, as illustrated in Figure 5a,b. The initial velocity was set to be in the Z direction, which is normal to the rigid block. Since the shell structure is not a closure surface, there should be a concentrated force applied together with the added mass to resist the external pressure. The total force value was 2 q R  . Figure 6a,b shows the details of the meshes for both models. There are 7935 shell elements and 2600 solid elements after optimal meshing. Since the pressure load exists at all times, the static analysis under the pressure will only be carried out as the precursor prior to the dynamic analysis in LS-DYNA [56][57][58]. The external pressure can be calculated by Equation (6): where ρ, g, and h are the density of water and gravity acceleration and the height under the sea, respectively. The pressure is then applied to the surface of the shell, as illustrated in Figure 5a,b. The initial velocity was set to be in the Z direction, which is normal to the rigid block. Since the shell structure is not a closure surface, there should be a concentrated force applied together with the added mass to resist the external pressure. The total force value was qπR 2 . Figure 6a,b shows the details of the meshes for both models. There are 7935 shell elements and 2600 solid elements after optimal meshing. Since the pressure load exists at all times, the static analysis under the pressure will only be carried out as the precursor prior to the dynamic analysis in LS-DYNA [56][57][58].  Figure 7 shows the initial state of X-stress, which is caused by the static pressure only without the time effect. Figure 7a shows the X-stress of the cylinder model, while Figure 7b shows that of the spherical shell model. It can be found that the maximum stress for each model is −648 MPa and −660 MPa (the negative sign implies that the member is in a compression state), respectively, which only has a 2.0% discrepancy. In addition, the X-stress at the apex will be    Figure 7 shows the initial state of X-stress, which is caused by the static pressure only without the time effect. Figure 7a shows the X-stress of the cylinder model, while Figure 7b shows that of the spherical shell model. It can be found that the maximum stress for each model is −648 MPa and −660 MPa (the negative sign implies that the member is in a compression state), respectively, which only has a 2.0% discrepancy. In addition, the X-stress at the apex will be qR/2t = 650 MPa. The result agrees very well.  Figure 7 shows the initial state of X-stress, which is caused by the static pressure only without the time effect. Figure 7a shows the X-stress of the cylinder model, while Figure 7b shows that of the spherical shell model. It can be found that the maximum stress for each model is −648 MPa and −660 MPa (the negative sign implies that the member is in a compression state), respectively, which only has a 2.0% discrepancy. In addition, the X-stress at the apex will be       Figure 10a depicts the time history of the kinetic energy. It shows that the kinetic energy is almost the same for both cases up until the half of the contact time (t < 0.045 s); then, there is some time delay (phase difference) between the cases of the cylinder model and the spherical model. After the impulse contact, the values of kinetic energy remain mostly the same (t > 0.08). Figure 10b depicts the internal energy and total energy time history of both models. It could be noticed that the internal energy of the cylinder is two times of that of the spherical shell model before collision (t < 0.02). Note that the cylinder is set to be the rigid material, so that it cannot absorb energy [59]. Nevertheless, the changing trends of the internal energy and total energy are the same.     Figure 10a depicts the time history of the kinetic energy. It shows that the kinetic energy is almost the same for both cases up until the half of the contact time (t < 0.045 s); then, there is some time delay (phase difference) between the cases of the cylinder model and the spherical model. After the impulse contact, the values of kinetic energy remain mostly the same (t > 0.08). Figure 10b depicts the internal energy and total energy time history of both models. It could be noticed that the internal energy of the cylinder is two times of that of the spherical shell model before collision (t < 0.02). Note that the cylinder is set to be the rigid material, so that it cannot absorb energy [59]. Nevertheless, the changing trends of the internal energy and total energy are the same.  Figure 10a depicts the time history of the kinetic energy. It shows that the kinetic energy is almost the same for both cases up until the half of the contact time (t < 0.045 s); then, there is some time delay (phase difference) between the cases of the cylinder model and the spherical model. After the impulse contact, the values of kinetic energy remain mostly the same (t > 0.08). Figure 10b depicts the internal energy and total energy time history of both models. It could be noticed that the internal energy of the cylinder is two times of that of the spherical shell model before collision (t < 0.02). Note that the cylinder is set to be the rigid material, so that it cannot absorb energy [59]. Nevertheless, the changing trends of the internal energy and total energy are the same.    Figure 10a depicts the time history of the kinetic energy. It shows that the kinetic energy is almost the same for both cases up until the half of the contact time (t < 0.045 s); then, there is some time delay (phase difference) between the cases of the cylinder model and the spherical model. After the impulse contact, the values of kinetic energy remain mostly the same (t > 0.08). Figure 10b depicts the internal energy and total energy time history of both models. It could be noticed that the internal energy of the cylinder is two times of that of the spherical shell model before collision (t < 0.02). Note that the cylinder is set to be the rigid material, so that it cannot absorb energy [59]. Nevertheless, the changing trends of the internal energy and total energy are the same. According to the discussion above, it can be seen that all of the physical quantities of force, stress, displacement, and energy of both models could be the same. Thus, it is necessary to consider the spherical model instead of a simplified cylinder model in order to determine the dynamic buckling of the structure subjected to the collision impacts.

Relationship between Static Force and Dynamic Force
From Section 2 above, the static buckling force, including point load and pressure load, can be calculated by Equation (6). To identify the buckling under the dynamic force of collision impacts, the relationship between the static force and dynamics force should be established.

Methodology
In this section, the relationship between the static and dynamic force is investigated. Figure 11 illustrates the schematic plot of the calculation for the impulse, which is referred to as Im , the average force, which is referred to as m f , and the max force, which is referred to as p f . The average force and impulse can be calculated by Equation (7) as follows: where 1 t and 2 t are the beginning and end of the collision procedure, According to the discussion above, it can be seen that all of the physical quantities of force, stress, displacement, and energy of both models could be the same. Thus, it is necessary to consider the spherical model instead of a simplified cylinder model in order to determine the dynamic buckling of the structure subjected to the collision impacts.

Relationship between Static Force and Dynamic Force
From Section 2 above, the static buckling force, including point load and pressure load, can be calculated by Equation (6). To identify the buckling under the dynamic force of collision impacts, the relationship between the static force and dynamics force should be established.

Methodology
In this section, the relationship between the static and dynamic force is investigated. Figure 11 illustrates the schematic plot of the calculation for the impulse, which is referred to as Im, the average force, which is referred to as f m , and the max force, which is referred to as f p . The average force and impulse can be calculated by Equation (7) as follows: where t 1 and t 2 are the beginning and end of the collision procedure, ∆t = t 2 − t 1 , and f(t) is the contact force, which is the variable with respect to time during the contact. The skewness of the contact force is measured by k s = f p / f m . The ratio between the max and mean force represents the dynamic contact load redistribution or the skewness of the dynamic collision force during the impact. According to the discussion above, it can be seen that all of the physical quantities of force, stress, displacement, and energy of both models could be the same. Thus, it is necessary to consider the spherical model instead of a simplified cylinder model in order to determine the dynamic buckling of the structure subjected to the collision impacts.

Relationship between Static Force and Dynamic Force
From Section 2 above, the static buckling force, including point load and pressure load, can be calculated by Equation (6). To identify the buckling under the dynamic force of collision impacts, the relationship between the static force and dynamics force should be established.

Methodology
In this section, the relationship between the static and dynamic force is investigated. Figure 11 illustrates the schematic plot of the calculation for the impulse, which is referred to as Im , the average force, which is referred to as m f , and the max force, which is referred to as p f . The average force and impulse can be calculated by Equation (7) as follows: The ratio between the max and mean force represents the dynamic contact load redistribution or the skewness of the dynamic collision force during the impact.  Table 4 presents the statistical data of the dynamic forces under different velocities. The results of ∆t, the maximum force, and the average force are calculated as described above.  Figure 12 shows the relationships between forces and velocity simulated under various cases of thickness at 0.04 m and 0.06 m. Figure 13 shows the time history of forces for the cases of impact velocities of 12 m/s, 14 m/s, and 16 m/s, respectively. It should be noted that when the velocity is larger than 14, the shell structure can be considered to be buckled already. According to Figure 13 and Table 4, the duration and impulse of collision contact is almost the same before buckling, and the skewness ratio of k s is almost the same as well, which implies that the contact forces have the same shape functions under a variety of impact velocities.

Parametric Effect Due to Velocity
Appl. Sci. 2018, 8, x 11 of 23 Figure 11. The typical plot of force vs. time. Table 4 presents the statistical data of the dynamic forces under different velocities. The results of t  , the maximum force, and the average force are calculated as described above.   Figure 12 shows the relationships between forces and velocity simulated under various cases of thickness at 0.04 m and 0.06 m. Figure 13 shows the time history of forces for the cases of impact velocities of 12 m/s, 14 m/s, and 16 m/s, respectively. It should be noted that when the velocity is larger than 14, the shell structure can be considered to be buckled already. According to Figure 13 and Table 4, the duration and impulse of collision contact is almost the same before buckling, and the skewness ratio of s k is almost the same as well, which implies that the contact forces have the same shape functions under a variety of impact velocities.    Table 4 presents the statistical data of the dynamic forces under different velocities. The results of t  , the maximum force, and the average force are calculated as described above.   Figure 12 shows the relationships between forces and velocity simulated under various cases of thickness at 0.04 m and 0.06 m. Figure 13 shows the time history of forces for the cases of impact velocities of 12 m/s, 14 m/s, and 16 m/s, respectively. It should be noted that when the velocity is larger than 14, the shell structure can be considered to be buckled already. According to Figure 13 and Table 4, the duration and impulse of collision contact is almost the same before buckling, and the skewness ratio of s k is almost the same as well, which implies that the contact forces have the same shape functions under a variety of impact velocities.   In addition, it can be observed that the maximum force and the mean force during collision are proportional to the velocity before dynamic buckling incurs, and the impulse is also proportional to the velocity. To identify the skewness ratio, both the maximum force and the mean force can be written in the same function as Equation (8) with different coefficients:

Parametric Effect Due to Velocity
where k v and km v are the proportional ratio. The results of regression analysis are shown in Table 5.
The skewness can be descripted by the ratio of k s = k v /km v , which can be considered as a constant value of 1.70. On this ground, it can be concluded that the maximum force and the average force comply with the same rules, and the ratio is dependent on the thickness. Moreover, the ratios of k v and km v in different cases of t = 0.06 m and 0.04 m are 1.5 and 1.6, respectively, which are almost the same as the identical ratio over the corresponding thickness. It is important to note that Equation (8) is independent of the other parameters.

Parametric Effect Due to Thickness
The  Table 6 presents the relationship between the contact force and the thickness. It should be noted that the impact duration is different in each case of thickness, since the shell thickness can change the contact stiffness. To compare the effect of the shell thickness, the contact forces due to an impact velocity of 1.0 m/s are visualized over the data derived from an impact velocity of 4.0 m/s, as illustrated by Equation (6) and Table 4 (previously discussed in Section 3.3.2).  Figure 14a shows the relationship between force versus impulse and the shell thickness considering the velocity of 4.0 m/s. Figure 14b shows the same for the velocity of 8.0 m/s. It can be seen that the impulse is almost the same, despite the differences in impact duration. It can be found that the relationship can be expressed as a linear relation as written in Equation (9), considering the results from the various conditions shown in Table 7, where, k t and km t are the coefficients of maximum force and mean force, and C t and Cm t are a constant. When the impulses are almost identical, it is noted that the higher the max force, the shorter the time duration.  According to the results shown in Table 7, it can be found that the ratio of  Table 8 presents the relationship between the contact forces and elastic modulus ratios together with external pressures. It should be noted that the duration time of impulse contact decreases with the incremental increase of elastic modulus (or contact stiffness), but the duration time remains the same over a variety of external pressure values.   According to the results shown in Table 7, it can be found that the ratio of k t /km t to C t /Cm t can be considered as a fixed value at 1.70, which is the average value of 1.72, 1.69, 1.67 and 1.72; thereby, the maximum force and the mean force can be expressed using an identical formulation, f p = k s · f m . Furthermore, it can be established that the contact forces during collisions have an identical shape function, despite the differences in shell thickness. The ratios of k t to km t under different cases of v = 4.0 and 8.0 m/s are 0.47 and 0.48, respectively, which are almost identical to the ratio of the corresponding velocity. Thus, it can be evident that Equation (9) is independent of the other parameters.

Effect of Elastic Modulus and External Pressure
Numerical simulations have been carried out for evaluating the effects of elastic modulus and external pressure. The shell models have adopted the radius of 2.0 m, a Poisson ratio of 0.15, the external pressure of 15 MPa, the velocity of 3.0 m/s, and thicknesses of 0.04 m and 0.05 m, with a variety of elastic moduli and different pressures, respectively. Note that e is the elastic modulus ratio, which is defined as e= E/210. Table 8 presents the relationship between the contact forces and elastic modulus ratios together with external pressures. It should be noted that the duration time of impulse contact decreases with the incremental increase of elastic modulus (or contact stiffness), but the duration time remains the same over a variety of external pressure values.  Figure 15a illustrates the relationships between force and elastic modulus ratios, while Figure 15b illustrates the relationships between forces and external pressures. Based on the results obtained, the contact forces comply with Equation (10), while the impulses remain the same over different cases. These results imply that the exchange of impact momentum for every case under various elastic modulus ratios or different external pressures is identical. f p = k e · e + C e , f m = km e · e + Cm e (10) where k e and km e are the proportional coefficients. The results of regression analysis are shown in Table 9. It is apparent that the ratio of k e /km e and C e /Cm e can be considered as a constant value of 1.70, and the slope of the force over the pressure relationship is relatively very small. It can be evident that the dynamic forces are constant, regardless of the changes in external pressure. The maximum force and the mean force can be expressed in accordance with the analytical formulation, f p = k s · f m , where k s is 1.70 without the effect of the elastic modulus ratio. It is also clear that the contact forces during collision impacts have an identical shape function for all of the cases of the various elastic modulus ratios.  Figure 15a illustrates the relationships between force and elastic modulus ratios, while Figure  15b illustrates the relationships between forces and external pressures. Based on the results obtained, the contact forces comply with Equation (10), while the impulses remain the same over different cases. These results imply that the exchange of impact momentum for every case under various elastic modulus ratios or different external pressures is identical.
f k e C f km e Cm (10) where e k and e km are the proportional coefficients. The results of regression analysis are shown in Table 9. It is apparent that the ratio of  Table 10 presents the relationship between forces and mass ratios. It should be noticed that the duration time of contact increases with the incremental increase of the mass ratio. From the regression analysis, it can be obtained that the maximum force and the mean force can be expressed in Equation

Parametric Effect of Added Mass
The simulations have been carried out using the shell radius of 2.0 m, a Poisson ratio of 0.15, the external pressure of 15 MPa, the velocity of 3.0 m/s, and thicknesses of 0.04 m and 0.06 m with different added mass ratios. It should be noticed that the mass ratio can be calculated as Equation (11), where M is the added mass, and M 0 = 35,800 kg. The different mass ratio can simulate the different weights of the submarine. m = M/M 0 (11) Table 10 presents the relationship between forces and mass ratios. It should be noticed that the duration time of contact increases with the incremental increase of the mass ratio. From the regression analysis, it can be obtained that the maximum force and the mean force can be expressed in Equation (12), where the k M , C M , km M , and Cm M values are the proportional coefficients respective to the mass ratios.  Figure 16a illustrates the relationship between the maximum force and the mean force over a variety of mass ratios, and Figure 16b illustrates the relationship of both forces versus external pressures. It is apparent that both the maximum force and the mean force comply with a linear function, as expressed by Equation (12). In addition, the ratios of k M /km M and C M /Cm M are almost the same value of 1.70. It can be observed that in the cases of different added masses, the contact forces have an identical skewness ratio. Furthermore, the ratios of k M and km M in different thickness cases of t = 0.04 meters and t = 0.06 meters are 0.6 and 0.66 respectively, which are almost identical to the ratios of the corresponding thicknesses. Thus, it can be concluded that Equation (12) is independent of the other parameters.  Figure 16a illustrates the relationship between the maximum force and the mean force over a variety of mass ratios, and Figure 16b illustrates the relationship of both forces versus external pressures. It is apparent that both the maximum force and the mean force comply with a linear function, as expressed by Equation (12). In addition, the ratios of  (12) is independent of the other parameters.

Contact Force Principle
Based on the results obtained earlier, the average force during the collision contact is proportional to the thickness and velocity, and it is linear to the elastic modulus and added mass based on the elastic material. Thus, the force can be written as Equation (13): Combining equations (8) and (12), and taking into account the data in Tables 4-11, the formula can be established precisely, and it can be used to solve for a variable when knowing the other parameters. Simultaneously, the collision contact forces have the almost same skewness ratio during collision. Thus, the mean force can be calculated by Equation (14)

Contact Force Principle
Based on the results obtained earlier, the average force during the collision contact is proportional to the thickness and velocity, and it is linear to the elastic modulus and added mass based on the elastic material. Thus, the force can be written as Equation (13): Combining equations (8) and (12), and taking into account the data in Tables 4-11, the formula can be established precisely, and it can be used to solve for a variable when knowing the other parameters. Simultaneously, the collision contact forces have the almost same skewness ratio during collision. Thus, the mean force can be calculated by Equation (14): In addition, since the contact forces in various cases using different impact velocities have the same shape function, the contact forces can thereby be descripted as Equation (15). Figure 17 shows the comparison of contact forces between the fitting results derived from Equation (15) and the numerical data. Excellent agreement can be established. Besides, according to Equation (15), the skewness ratio can be derived theoretically as k s = π/2 = 1.57, which has a 7.6% discrepancy compared with the numerical value 1.70.
In addition, since the contact forces in various cases using different impact velocities have the same shape function, the contact forces can thereby be descripted as Equation (15). Figure 17 shows the comparison of contact forces between the fitting results derived from Equation (15) and the numerical data. Excellent agreement can be established. Besides, according to Equation (15), the skewness ratio can be derived theoretically as /2 s k  = = 1.57, which has a 7.6% discrepancy compared with the numerical value 1.70. Table 11. The result of regression analysis on Equation (12). Ks f t t (16) It can be evident that the function of displacement versus time is a harmonic function as well, and the duration time of collision depends on the thickness, elastic modulus, and density of the shell, but has little effect due to the radius and velocity of the spherical shell. This conclusion can be verified by the extensive data in tables 4-10.

Dynamic Buckling Criteria
The dynamic buckling problem can be solved if the relationship between the static buckling force and maximum force is established. In this section, the relationships between the maximum forces of collision and static buckling forces according to Equation (3) are investigated. Figure 18 shows the X-stress of the shell under the external pressure of 15 MPa. At the beginning, the initial state is the same between the cases of velocities of 2.8 m/s and 2.7 m/s when the pressure load and configuration of the shell are the same. According to Equation (15), the dynamics equation can be written in Equation (16) correspondingly, where K is the stiffness of the spherical shell. The relationship depends on the thickness, elastic modulus, and radius of the shell: It can be evident that the function of displacement versus time is a harmonic function as well, and the duration time of collision depends on the thickness, elastic modulus, and density of the shell, but has little effect due to the radius and velocity of the spherical shell. This conclusion can be verified by the extensive data in Table 4, Table 5, Table 6, Table 7, Table 8, Table 9, Table 10.

Dynamic Buckling Criteria
The dynamic buckling problem can be solved if the relationship between the static buckling force and maximum force is established. In this section, the relationships between the maximum forces of collision and static buckling forces according to Equation (3) are investigated. Figure 18 shows the X-stress of the shell under the external pressure of 15 MPa. At the beginning, the initial state is the same between the cases of velocities of 2.8 m/s and 2.7 m/s when the pressure load and configuration of the shell are the same. Figure 19a-h shows the mechanisms of collision for the cases with impact velocities of 2.8 m/s and 2.7 m/s. It demonstrates that in the velocity case of 2.7 m/s, the deformation of the spherical shell can be restored elastically after the collisions. However, when the velocity is larger than 2.7 m/s (even at 2.8 m/s), the deformation of the spherical shell cannot be restored after the collisions, which became larger and larger under the load of pressure and contact force. In that case, the shell can be considered to have buckled under the dynamic collisions.    Figure 19a-h shows the mechanisms of collision for the cases with impact velocities of 2.8 m/s and 2.7 m/s. It demonstrates that in the velocity case of 2.7 m/s, the deformation of the spherical shell can be restored elastically after the collisions. However, when the velocity is larger than 2.7 m/s (even at 2.8 m/s), the deformation of the spherical shell cannot be restored after the collisions, which became larger and larger under the load of pressure and contact force. In that case, the shell can be considered to have buckled under the dynamic collisions.  Figure 20 illustrates the Z displacement of three nodes, as previously shown in Figure 19h, which are marked as yellow spots with velocities of 2.8 m/s and 2.7 m/s, respectively. It can be found that the displacements are the same before collision contact, but they can be very different from the beginning of contact. In the velocity case of 2.7 m/s, the displacement will return back symmetrically, whilst in the case of 2.8 m/s, the displacement will return back in a steep slope (non-elastic), which implies that the shell had buckled already (the same conclusion as with Figure 19).

Dynamic Buckling Force
On this ground, the critical impact velocity Vcr can be determined to be 2.7 m/s, and according to the criteria shown above, some more critical velocities in different cases can be obtained in Table  Figure 20 illustrates the Z displacement of three nodes, as previously shown in Figure 19h, which are marked as yellow spots with velocities of 2.8 m/s and 2.7 m/s, respectively. It can be found that the displacements are the same before collision contact, but they can be very different from the beginning of contact. In the velocity case of 2.7 m/s, the displacement will return back symmetrically, whilst in the case of 2.8 m/s, the displacement will return back in a steep slope (non-elastic), which implies that the shell had buckled already (the same conclusion as with Figure 19).  Figure 20 illustrates the Z displacement of three nodes, as previously shown in Figure 19h, which are marked as yellow spots with velocities of 2.8 m/s and 2.7 m/s, respectively. It can be found that the displacements are the same before collision contact, but they can be very different from the beginning of contact. In the velocity case of 2.7 m/s, the displacement will return back symmetrically, whilst in the case of 2.8 m/s, the displacement will return back in a steep slope (non-elastic), which implies that the shell had buckled already (the same conclusion as with Figure 19).

Dynamic Buckling Force
On this ground, the critical impact velocity Vcr can be determined to be 2.7 m/s, and according to the criteria shown above, some more critical velocities in different cases can be obtained in Table   Figure

Dynamic Buckling Force
On this ground, the critical impact velocity V cr can be determined to be 2.7 m/s, and according to the criteria shown above, some more critical velocities in different cases can be obtained in Table 12. It should be noted that for the application of submarines, an exploration is usually among 1.0 km to 4.5 km deep, and the external pressure is correspondingly set to be within 10 MPa to 45 MPa. According to Table 12, the critical load in the collision can be converted into the normalized results tabulated in Table 13. Note that the labels p' and q' show the normalized values of pressure and maximum force respectively, according to Equation (6). The p e shows the static buckling force calculated from Equation (6) as well. The ratio is deducted by p/p', which demonstrates the multiple of the maximum force over to the static buckling force from Equation (6). It can be found that the maximum critical force is about 2.1 times that of the static force, inducing the dynamic buckling during the collision impacts. In summary, the critical dynamic force over the static buckling force ratios can be illustrated in Figure 21. It can also be concluded that the analytical prediction fits with the numerical simulation result very well. 12. It should be noted that for the application of submarines, an exploration is usually among 1.0 km to 4.5 km deep, and the external pressure is correspondingly set to be within 10 MPa to 45 MPa.  It can be found that the maximum critical force is about 2.1 times that of the static force, inducing the dynamic buckling during the collision impacts. In summary, the critical dynamic force over the static buckling force ratios can be illustrated in Figure 21. It can also be concluded that the analytical prediction fits with the numerical simulation result very well.

Discussion
According to Equation (5), when the pressure tends to be zero, which means that there is no external pressure under the shell, the buckling forces can increase without boundary. Physically speaking, if the pressure ratio is less than 0.05, the force ratio would be larger than 2.8. In such a case, it can be considered that the structural buckling of the shell cannot occur.
During the collision impacts, the maximum force can be considered to be 1.70 times the average force derived from Equation (7). The impulse of collision is linear to the velocity and added mass, but has a negligible change due to the thickness and elastic modulus. In addition, the quasi-static buckling force deduced by Equation (5) is about only 0.48 times that of the maximum force during the collision. It should be noted that Equation (5) is based on the elastic material, as the buckling usually occurs before the plastic stage. The effect of strain rate is neglected for the low velocity

Discussion
According to Equation (5), when the pressure tends to be zero, which means that there is no external pressure under the shell, the buckling forces can increase without boundary. Physically speaking, if the pressure ratio is less than 0.05, the force ratio would be larger than 2.8. In such a case, it can be considered that the structural buckling of the shell cannot occur.
During the collision impacts, the maximum force can be considered to be 1.70 times the average force derived from Equation (7). The impulse of collision is linear to the velocity and added mass, but has a negligible change due to the thickness and elastic modulus. In addition, the quasi-static buckling force deduced by Equation (5) is about only 0.48 times that of the maximum force during the collision. It should be noted that Equation (5) is based on the elastic material, as the buckling usually occurs before the plastic stage. The effect of strain rate is neglected for the low velocity analysis in this study. It is noted that the spherical shell could reach up to the plastic stage or over the limited stress under the collision before buckling. In that case, the structure will be damaged or failed for the lack of strength, so it is suggested that the static analysis should be analyzed before the buckling analysis.
For applications from the outcome of this study, the critical force under a single point load and a single uniform pressure load, which are named p c and q c , respectively, can be calculated firstly by Equation (1) and Equation (2) relevant to the property of the material and the parameters of the geometry. Then, the external pressure q can be determined by the operation depth under the sea; therefore, Equation (5) can be used by the designer to determine the quasi-static buckling force. After obtaining the normalized force ratio of p , it can be used to identify the critical velocity, which may induce the dynamic buckling by Equation (13). By the mechanism described above, the dynamic buckling of the spherical shell under the collision impacts and external pressure can be evaluated.

Conclusions
This paper has studied the buckling of a spherical shell structure subjected to both external pressure and point load in order to establish criteria to predict the dynamic buckling of a spherical shell structure under collision impacts. In the present study, the standard for static buckling, including external pressure and concentrated load, was established as Equation (5). Furthermore, the relationship between static force and dynamic force under subsea collisions has been studied based on the theory that the static buckling criteria can be used to reflect the dynamic collision problem. By extensive parametric analyses, a precise formulation has been formed and validated using multiple regression analyses. During the collisions, the maximum and mean contact forces are proportional to the velocity and thickness of the shell, and are also linear to the total mass and elastic modulus. For the collision impulses, they are proportional to the velocity and added mass, but the elastic modulus and shell thickness do not affect their momentum change. In addition, the same skewness ratio of the contact forces during different parameters can be considered the same value: 1.70. In that case, the contact force during collision can be descripted as a sine function with respect to time.
Overall, this study has established a fundamental framework and novel theory to simulate, design, and validate the dynamic buckling for a spherical shell structure under collision impacts and external pressure. The precise formulations shown in the present paper enable the prediction of the dynamic buckling of spherical shell structures subjected to subsea collisions. Such insights and outcomes will benefit the potential applications involving submarines and deep sea exploration, etc.
Author Contributions: L.P., S.K. and Z.D.C. developed and initialized the concept and methods; L.P. developed model and carried out analyses; S.K. provided analysis tools and validated the model; Z.D.C. and S.S.W. provided industry guidance and verified results with national standards; L.P., S.K., Z.D.C. and S.S.W. wrote and reviewed the paper; L.P. and S.S.W sought funding.
Funding: This work was supported by the national Natural Science Foundation of China (grant numbers 51508238) and by the Jiangsu Postdoctoral Research Plan (grant numbers 1601014B). The APC is sponsored by the University of Birmingham Library, U.K.