A meshfree framework for the numerical simulation of elasto-plasticity deformation of ship structure

Ocean Engineering 192 (2019) 106507

Contents lists available at ScienceDirect

Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng

A meshfree framework for the numerical simulation of elasto-plasticity deformation of ship structure Yu-Xiang Peng, A-Man Zhang *, Fu-Ren Ming, Shi-Ping Wang College of Shipbuilding Engineering, Harbin Engineering University, Harbin, 150001, China

A R T I C L E I N F O

A B S T R A C T

Keywords: Elasto-plasticity Meshfree methods Ship hull Dynamic response Large deformation

In present work, a meshfree framework for the numerical simulation of elasto-plasticity deformation of ship structure is established. The reproducing kernel particle method (RKPM) is adopted to discretize the calculation domain, and the meshfree shell formulation is derived based on Mindlin shell theory. Then a general elastoplasticity model is established under the meshfree shell framework, and the stress resultants are avoided since the shell formulation is based on degenerated continuum approach. By removing the normal stress from the yield function, the zero normal stress assumption is implemented. Several numerical examples are studied and the accuracy of the present elasto-plasticity model is verified. Finally, the dynamic response of the ship cabin and real ship structure under impact load are calculated, the capability of the present meshfree approach applied in complex engineering structures analysis is demonstrated.

1. Introduction Shell structures have been widely applied in many aspects of engi­ neering, especially in shipbuilding engineering, since the main structure of the ship hull is built by stiffened shell structures. Therefore, the response of the stiffened shell structure under ultimate load is of great concern to the engineers, for instance, the car crash (Cheng et al., 2001), ship grounding (Simonsen, 1997), and damage of ship hull subjected to the underwater explosion (Liang and Tai, 2006; Ming et al., 2016). All these cases contain large deformation of the shell structures. In the numerical simulation of such problems, the traditional finite element method (FEM) is facing difficulties, such as mesh distortion and dynamic crack propagation. However, the meshfree methods have an advantage in handling these cases, and they do not suffer from mesh distortion since their discretizations do not rely on the finite mesh (Li and Liu, 2002; Nguyen et al., 2008). In recent studies, the meshfree methods have been applied in the dynamic fracture simulation of thin shells. For instance, Rabczuk et al. (2007) proposed a meshfree thin shell method for the crack propagation in thin shells; Rabczuk et al. (2010) using the meshfree method dealing with dynamic fracture of shells while the fluid-structure interaction is considered. For the numerical simulation of shell structures, there are two main

approaches: one is the direct three-dimensional continuum approach (Li et al., 2000); the other one is the so-called degenerated continuum approach (Wang and Chen, 2004; Rabczuk and Belytschko, 2004; Areias and Rabczuk, 2013). The numerical implementation of the former approach is simpler than the degenerated continuum approach and al­ lows an implementation of local three-dimensional constitutive model without any modifications. However, the former approach is not popular in practice, since it has difficulties to model the complex shell structures and the computational cost is relatively high. In contrast, the degen­ erated continuum approach is not limited by these problems and has been widely used in numerical simulation (Kiendl et al., 2009; Amiri et al., 2014; Nguyen-Thanh et al., 2015, 2017). Moreover, the deep neural network (DNN) is applied in the bending analysis of thin shells recently, related work on this reader may refer to Guo et al. (2019). In previous research, the meshfree methods have been applied in the simulation of shell structures. Caleyron et al. (2012) applied the Moving least squares particle hydrodynamics (MLSPH) to the simulation of the dynamic fracture of shell structures. Ren and Li (2012) used the repro­ ducing kernel particle method (RKPM) in the simulation of large-scale ductile fracture of plates and shells; O’Grady and Foster (2014) devel­ oped a peridynamics shell formulation to model the bending deforma­ tion. However, the meshfree methods have seldom been applied to the analysis of complex shell structure, for instance, the meshfree methods

* Corresponding author. E-mail addresses: [email protected] (Y.-X. Peng), [email protected] (A.-M. Zhang), [email protected] (F.-R. Ming), shipingwang316@126. com (S.-P. Wang). https://doi.org/10.1016/j.oceaneng.2019.106507 Received 20 May 2019; Received in revised form 19 September 2019; Accepted 28 September 2019 Available online 16 October 2019 0029-8018/© 2019 Elsevier Ltd. All rights reserved.

Y.-X. Peng et al.

Ocean Engineering 192 (2019) 106507

have not been applied to the simulation of a real ship structures. Focus on the insufficiency of current research, in the present work, we estab­ lished a meshfree framework for the numerical simulation of elasto­ plasticity response of ship structures. In present work, the RKPM shell formulation based on Mindlin shell theory Peng et al. (2018) are adopted, which characterized by the straight fibers remain within the deformation process. The elastoplastic constitutive model applied in degenerated shell has been studied by many researchers (Simo and Kennedy, 1992; Areias and Belytschko, 2006; Klinkel et al., 2008). The model proposed by Simo and Kennedy (1992) is efficient and has been simplified by Zeng et al. (2001) to further improve the computational efficiency. However, plastic model based on stress resultants is less hardening than the classical J2 model. In this work, we draw on the matrix analysis method which has been adopted by Simo and Kennedy (1992). But the use of stress resultants is avoided, the stress tensor of the integration point is chosen as the input variable of the plastic model instead. The yield criterion is derived from the von Mises yield criterion, and it can easily combine with the damage model further. In addition, since the stress resultants is not adopted, the dimension of the matrix used in the return mapping algorithm has reduced from 8 to 5, thus the complexity of the numerical algorithms is decreased. Another problem to solve is the implementation of the so-called zero normal stress condition, which is essential for the simulation of shell structures. In the literature, at the end of each iteration, the thickness is changed to satify the zero normal stress condition. In this work, by considering this condition in the yield criterion and using the reduced constitutive model, the zero normal stress condition is satisfied. The shell simulation is based on the RKPM shell formulation proposed by Peng et al. (2018), which has been used for the dynamic analysis of stiffened shell structures (Peng et al., 2019). The paper is organized as follows. A brief introduction of the RKPM shell formulation and the reproducing kernel function is presented in section 2. In section 3, the yield criterion combined with the zero normal stress assumption is present, and the return mapping algorithm is derived. In section 4 we present several numerical examples to verify the reliability of the present model. Lastly, a summary of this paper is pre­ sented in section 5.

vðξ; η; ζÞ ¼

� X� h~ NI vI þ N I ωI � yI 2 I

(2)

~ I ¼ ζNI , vI is the nodal velocity vector, ωI is the fiber rotation where N angular velocity. The local corotational coordinate system at each stress point, feli ; i ¼ 1; 2; 3g, is defined by three mutually orthogonal unit vectors: x;ξ � el1 ¼ � �x;ξ �

(3)

el � x;η � el3 ¼ �� 1l je1 � x;η j�

(4)

el2 ¼ el3 � el1

(5)

where el1 ; el2 are the unit tangent vector of the reference surface, el3 is the normal of the reference surface. The local coordinate basis vector rotates rigidly while the shell deforms. The rotation matrix q that transform quantities from the global coordinate system to the local coordinate T

system is defined as q ¼ ½el1 ; el2 ; el3 � , where the superscript’ T’ denotes transpose. 2.2. RKPM shape function In Eq. (1), the interpolation shape function NI ðξÞ is introduced. In present work, the reproducing kernel function is adopted, which is constructed by correcting the original SPH kernel function. Besides, the shape function only calculated once at the begin of the simulation, since the Lagrangian kernel is adopted. The shape function used in RKPM shell formulation can be written as (Peng et al., 2018): NI ðξÞ ¼ Cðξ; ξI

ξÞφa ðξI

ξÞΔξI

(6)

where ξI ξ ¼ q0 ⋅ðXI XÞ is the relative position between two neighbor points in the laminar coordinate system, where q0 is the initial rotation matrix. φa ðξI ξÞ is the so-called window function, ΔξI denotes the area of node I on the reference surface, and Cðξ; ξI ξÞ is the correction function:

2. Outline of RKPM shell formulation

Cðξ; ξI

2.1. Shell kinematics

ξÞ ¼ pT ðξI

ξÞbða; ξÞ

(7)

where pT denotes the interpolation basis function, and the bilinear basis function is adopted. bða; ξÞ is a coefficient vector, which is determined by the reproducing condition (Chen et al., 1996), we have:

Based on the degenerated solid element approach developed by Hughes and Liu (1981), the kinematics of the shell are defined, as shown in Fig. 1. Therefore, the geometric definition of the shell is reduced. Specifically, a reference surface and the fiber vectors vertical to it are chosen to describe the shell kinematics. In numerical practice, only one layer of nodes on the reference surface is needed to discretize the shell. For the Galerkin type weak form, the field variable is reproduced by values of nodes in a generally domain Ω0 , and the shape function is used to approximate the field variable. For instance, in the current configu­ ration, we can express the position vector of a generic point as follows: � X � h xðξ; η; ζÞ ¼ NI xI þ ζyI (1) 2 I

bða; ξÞ ¼ M 1 ða; ξÞpT ð0Þ

where NI denotes the interpolation shape function, xðξÞ is the current coordinate vector of a generic point in the shell, xI is the current position vector of node I on the reference surface, and yI is the current nodal unit fiber vector, which is not required to be normal to the reference surface since the Mindlin-Reissner shell theory is adopted. Moreover, a rotation operator based on the quaternion is used to update the fiber vector. Thus, the velocity of a generic point in the shell can be obtained from the time derivative of the position vector:

(8)

Fig. 1. Kinematic description and the local coordinate system of RKPM shell Peng et al. (2018). 2

Y.-X. Peng et al.

Ocean Engineering 192 (2019) 106507

Z Mða; ξÞ ¼

ξÞpT ðξ’

pðξ’

ξÞφa ðξ’

Γ0

ξÞdξ’

Δt:

(9)

uptþΔt ¼ ut þ Δtu_ t þ

where pT ð0Þ ¼ ½1; 0; 0; 0�, pT ðξ’ ξÞ ¼ ½1;ξ’ ξ; η’ η;ðξ’ ξÞðη’ ηÞ�, and ξ’ dentoes and arbitrary point in the support domain. When the shell is discretized into nodes, the integration of Eq. (9) is calculated by the summation, and the differential area dξ’ is changed to ΔξI .

u_ ptþΔt ¼ u_ t þ ð1

The linear momentum equation for a continuum in the reference configuration Ω0 is:

ρ0 bðX; tÞ

r0 ⋅ P ¼ 0

(10)

Fint I

(11)

_ J ¼ Mext JIJ ω I

Mint I

(12)

¼ Ω0

ρ0 NI bdΩ0 þ

Γt0

NI t

dΓ t0

qθ ¼ cos

(16)

Ω0

Z

ρ0 h ~ 2

N I yI � bdΩ0 þ

h~ 0 N I yI � t dΓt0

Γ t0 2

(18)

�

�θ� θ y sin j 2 θ

�θ� θ z sin k 2 θ

(27)

3. Elastoplastic model of shell In present work, the nonlinear plasticity model for shell developed by Simo and Kennedy (1992) has been modified. The stress tensor of the Gauss point in the shell thickness direction instead of the stress re­ sultants is chosen as the input of the constitutive algorithm.

To update the displacement and rotations, firstly, the linear and angular accelerations can be calculated by Eqs. (11) and (12). Then the Newmark time integration scheme is adopted to update the displace­ ment and rotations. For instance, in the linear momentum equation Eq. (11), the mass matrix has been diagonalized by the row-sum lumped technique, thus, the linear acceleration of node I can be calculated by: Fint I

�θ � θ x sin i 2 θ

(26)

where y1 , y2 and y3 denote the three components of the fiber vector y, θx , θy and θz denote the rotation angles around three axis, and θ ¼ jjθjj, i, j, k is the unit basis vector of the global coordinate system. It should be noticed that the product of two quaternions in Eq. (24) is different from both the dot and cross product of two vectors, for detail, readers may refer to a comprehensive introduction by Kuipers (1999).

2.4. Time integration scheme and essential boundary conditions

1 €I ¼ u Fext mI I

(25)

�θ� �θ � θ �θ � θ �θ � θ x y z þ sin i þ sin j þ sin k 2 2 θ 2 θ 2 θ

�θ � qθ 1 ¼ cos 2

where dΩ0 ¼ hdΓ 0 , and t denotes the traction force. And the internal and external moment vector in Eq. (12) are expressed as follows: Z ~I h ∂N Mint (17) P � yI dΩ0 I ¼ Ω0 2 ∂X Z

(24)

pt ¼ 0 þ y1 i þ y2 j þ y3 k

0

Mext I ¼

(23)

where p is a pure quaternion, which is related to the fiber vector, qθ is related to the rotation angle and is a unit quaternion, therefore, its in­ verse qθ 1 can be obtained simply by changing the sign of its imaginary components. We have:

The internal force and external force vector in Eq. (11) are given as: Z ∂NI Fint (15) PdΩ0 I ¼ Ω0 ∂X 0

u_ tþΔt ¼ u_ ptþΔt þ γΔt€ utþΔt

pt ¼ qθ p0 qθ 1

where Γ0 denotes the initial reference surface. Similarly, we can obtain the lumped rotational inertia matrix, which can be expressed as: Z h2 h2 JIJ ¼ ρ hNI δIJ dΓ 0 ¼ mIJ (14) 12 Γ0 0 12

Z

(22)

For the simulation of stiffened shells, a support domain cutting scheme is developed. In this scheme, the support domain of a stress point, not only influenced by the distance between the two points but also related to which component the two points located. For detail, reader may refer to Peng et al. (2019). In present work, the shell formulation is based on the Mindlin shell theory, therefore, the fiber direction should be updated explicitly. The rotation operator based on the quaternion is adopted, which can be expressed as follows:

Γ0

Z

€tþΔt utþΔt ¼ uptþΔt þ βΔt2 u

2.5. Simulation of stiffened shells and fiber update scheme

where m and J is the mass and rotational inertia matrix, F, M denote the force and moment vector respectively, in which the superscript’ ext’ and’ int’ is used to identify the external or internal variables. A lumped mass matrix is obtained based on the row-summing technique, we have: Z mIJ ¼ ρ0 hNI δIJ dΓ0 (13)

Fext I

(21)

γÞΔt€ ut

The same scheme is employed to update the angular displacements and velocities. For the enforcement of the essential boundary conditions, the D’Alembert’s principle Günther and Liu (1998) is adopted.

where P is the first Piola-Kirchhoff stress, ρ0 denotes the initial density, b denotes the body force, and vðx; tÞ denotes the velocity field. Based on the principle of virtual power, the weak form discrete governing equa­ tions (Peng et al., 2018) are derived: €J ¼ Fext mIJ u I

(20)

2βÞ€ ut

where Δt denotes the time increment, β and γ are parameters. In all simulations, β ¼ 0 and γ ¼ 0:5 are adopted. The strain and stress tensor can be calculated by the predicted deform configuration, and the in­ ternal force at t þ Δt can be obtained. Therefore, the acceleration at t þ Δt is got. Finally, the velocity and displacement at t þ Δt is corrected by:

2.3. Governing equation of RKPM shell

∂vðX; tÞ ρ0 ∂t

Δt2 ð1 2

3.1. Basic equations of the multi-surface elastoplastic model Consider that a node inital located at X moves to x in the current configuration, then the deformation gradient tensor is calculated by: � � ∂xi X ∂NI h ∂ζNI Fij ¼ (28) ¼ xIi þ yIi 2 ∂ Xj ∂Xj ∂Xj I

(19)

After the acceleration at t is obtained, the following formula (Peng et al., 2019) is employed to predict the velocity and displacement at t þ 3

Y.-X. Peng et al.

Ocean Engineering 192 (2019) 106507

projected onto the shell reference surface. For instance, the stress tensor has to be transformed into the local coordinate system σ l ¼ qT σ g q. In the sequel, all the calculations are carried out in the local coordinate system, for writing simplicity, superscript ’l’ is omitted. In present work, the standard von Mises yield criterion is adopted, which is suitable for many applications in metal plasticity. And the yield function can be expressed as: pffiffiffiffiffiffi φ ¼ 3J2 σ y (39)

where i; j ¼ 1; 2; 3. The determinant of F is denoted as J ¼ detðFÞ. The Almansi strain is calculated by the deformation gradient:

ε¼

1 I 2

F TF

1

�

(29)

In present work, the strain tensor is divided into elastic and plastic parts: ε ¼ εe þ εp based on an additive decomposition (Karaoulanis, 2013), where εe and εp denote the elastic and plastic part of the strain tensor respectively. The constitutive model for the so-called hyperelastic material is the simplest, in which the stress response is governed by a strain energy function W which depending on the elastic deformations. Therefore, the Helmholtz free energy potential ψ can be decomposed into the strain energy function Wðεe Þ and the hardening energy function which related to the internal variable (Eberlein and Wriggers, 1999):

ψ ðε; εp ; r; αÞ ¼ Wðεe Þ þ H ðr; αÞ

where J2 is the second deviatoric stress invariant, σ y denotes the yield stress and σy ¼ κðrÞ ¼ κ0 q, where κ0 denotes the initial yield stress. To eliminate the square root, Eq. (39) is rewritten as: φ ¼ ð3J2 κ2 Þ=κ0 . By considering the zero normal stress condition σ 33 ¼ 0, we have:

(30)

φ¼

where r and α denote the isotropic and kinematic hardening state vari­ able, respectively. Since the damage of material has not been taken into account, we have r ¼ p, where p denotes the accumulate plastic strain. Following the constraint of isotropy, the local dissipation inequality is obtained based on the second law of thermodynamics, in the absence of thermal terms, the dissipation function can be expressed: D ¼ σ : ε_ ψ_ � 0, and which can be expanded as: � � ∂W ∂H ∂H D¼ σ : α_ � 0 (31) : ε_ e þ σ : ε_ p ⋅ r_ ∂εe ∂r ∂α

∂H q¼ ; ∂r

∂H β¼ ∂α

φ ¼ σ T Aσ

D ðσ ; q; βÞ þ γ_ φðσ ; q; βÞ

ε_ p ¼ γ_

∂φ ∂σ

A¼

� 1 P κ0 0

(32)

1 W ¼ ðε 2

3 1 07 2 7 7 7 1 07 5 0 3

(43)

κ2 κ0

(44)

εp ÞT Cðε

εp Þ

(45)

2

W where C ¼ ∂ε∂ e ∂ε is the elastic tensor. Since the strain tensor is rewritten e

into a vector, the elastic tensor is reduced from a fourth-order tensor to a second-order tensor. Thus we have: 2 3 0 1 ν � � 6 7 E 6ν 1 Cn 0 0 7 C¼ ; Cn ¼ (46) 6 7 2 0 G s I2 1 ν 4 1 ν5 0 0 2

(36)

(38)

(42)

To express the stress strain relation into matrix form, similar to Eq. (42), the strain tensor is rewritten into column vectors: εT ¼ ½ε11 ; ε22 ; ε12 ; ε13 ; ε23 �. Then the strain energy function is assumed to have the following form:

(35)

∂φ α_ ¼ γ_ ∂β

6 1 � 6 0 6 ; P¼6 1 3I2 6 4 2 0

φ ¼ ðσ þ βÞT Aðσ þ βÞ

(34)

(37)

(40)

where In denotes the rank-n identity matrix. If the kinematic hardening is taken into account, then the back stress should be considered:

(33)

∂φ ∂q

r_ ¼ γ_

κ2 κ0

(41)

2

where γ_ is the so-called plastic multiplier. According to the Kuhn-Tucker theorem, maximum plastic dissipation can be obtained if γ_ satisfies: γ_ > 0; φ � 0; γ_ φ � 0

κ2 κ0

σ T ¼ ½σ11 ; σ22 ; σ12 ; σ13 ; σ23 �

By assuming the principle of maximum plastic dissipation, we consider a maximization problem with inequality constraints: maximize D subject to φðσ ; q; βÞ � 0, where φ is the yield function of plastic cri­ terion. This maximization problem under constraints can be solved by constructing the following Lagrange function: L ðσ ; q; β; γÞ ¼

��

where σ has been rewritten into a column vector, and A is a coefficient matrix, we have:

where q denotes the isotropic hardening stress variable, β denotes the uniaxial kinematic back stress tensor. Then the Clausius-Planck consti­ tutive inequality Eq. (31) takes the reduced form: D ¼ σ : ε_ p þ q⋅r_ þ β : α_ � 0

�

σ 11 ⋅ σ 22 þ 3σ213 þ 3σ 223

The equation above can be expressed into matrix form:

Then the constitutive relations for Cauchy stress and hardening variable are derived from Eq. (31) (Truesdell and Noll, 2004):

∂W σ ¼ e; ∂ε

1� 2 σ 11 þ σ222 þ 3σ212 κ0

where E is the Young’s modulus, and ν is the Poisson’s ratio. Gs ¼ Gks , in which G and ks are denote the shear modulus and shear correction factor respectively. Similarly, the hardening potential is also assumed to be quadratic:

where Eq. (35) is the loading/unloading conditions, and Eq. (36–38) is the evolution equations for the rate of plastic strain, isotropic and ki­ nematic hardening variables.

1 1 H ¼ Ep r2 þ αT Dα 2 2

(47)

where Ep and D ¼ 23 H’I5 denote the isotropic and kinematic hardening modulus, respectively. Substituting Eqs. (45) and (47) into Eq. (32), then we have:

3.2. Yield function under zero normal stress condition Based on the general framework above, we now establish the elas­ toplastic model for shell structure. Since the degenerated solid approach is adopted, the variables in the constitutive equations need to be 4

Ocean Engineering 192 (2019) 106507

Y.-X. Peng et al.

Fig. 2. Simply supported plate loaded by pressure impulse: problem description.

σ¼

∂W ¼ C⋅εe ; ∂εe

q¼

Ep r;

β¼

Dα

(48)

3.3. Elastoplastic return mapping algorithm

Fig. 4. Comparison FEM (ABAQUS).

The main algorithmic problem to be solved in this section is to up­ date the field fεn ; εpn ; rn ; αn g at tn to fεnþ1 ; εpnþ1 ; rnþ1 ; αnþ1 g at tnþ1 . Firstly,

βtrial nþ1 ¼

Cεpn

Dαn

κ2 ðrÞ ¼0 κ0

cost

between

RKPM

and

σ nþ1 ¼ σ trial nþ1

CΔεpnþ1

(52)

(49)

βnþ1 ¼ βtrial nþ1

DΔαnþ1

(53)

(50)

φnþ1 ðσ ; q; βÞ ¼ 0

(54)

where Δεpnþ1 and Δαnþ1 denote the increment of plastic strain and hardening variables in the n þ 1 time step, respectively. And they can be calculated by substituting Eq. (44) into Eqs. (36) and (38). Furthermore, we can observe that there is only one unknown variable γ nþ1 in Eq. (52–54). All of the following efforts are for solving γnþ1 , and the Newton iteration scheme is adopted. The yield function for the k-th iteration step is expressed as follows: ðkÞ � κ2 rnþ1 ðkÞ ðkÞ ðkÞ �T ðkÞ ðkÞ � φnþ1 ¼ σ nþ1 þ βnþ1 A σ nþ1 þ βnþ1 ¼0 (55) κ0

Substituting Eq. (49) into Eq. (44), we have: � � trial T trial trial trial φtrial nþ1 ¼ σ nþ1 þ βnþ1 A σ nþ1 þ βnþ1

the computational

yield surface, the nonlinear system of equations below has to be solved:

we assumed that εpnþ1 ¼ εpn and αnþ1 ¼ αn , then the trail stress and trial back stress at the n þ 1 time step are obtained:

σ trial nþ1 ¼ Cεnþ1

of

(51)

If φtrial nþ1 � 0, which indicates that the material is in an elastic loading or unloading state and there is no new plastic deformation occurs. Then trial we set ð�Þnþ1 ¼ ð�Þtrial nþ1 . If the φnþ1 > 0, which indicates that the material is in the plastic loading state. Thus, to make the stress state return to the

where the stress and back stress in k iteration step, and can be expressed as follows:

Fig. 3. Displacement histories of the center point and convergence to the reference solution.

Fig. 5. Roof subjected to the velocity impulse: problem description. 5

Y.-X. Peng et al.

Ocean Engineering 192 (2019) 106507

Fig. 8. Simply supported plate under uniform distributed load.

Fig. 6. Displacement history of gauging point.

�

�

�

4 3

��

ðkÞ k σ ðkÞ nþ1 þ βnþ1 ¼ I5 þ γ nþ1 2CA þ H’A

1 trial σ trial nþ1 þ βnþ1

�

(56)

where: �

� 2� 3 4 � 0 4 1 4 2Cn P þ H’P 5 3 2CA þ H’A ¼ 3 κ0 0 3Gs I2

(57)

The matrix C and P have the same characteristic subspaces (Simo and Kennedy, 1992). That is P ¼ QΛP QT ; Cn ¼ QΛCn QT where Q ΛCn is:

1

(58)

¼ QT is a orthogonal matrix, the diagonal matrices Λp and

Fig. 9. Load level versus deflection at the center of the plate.

Fig. 7. Deformed configuration compared to experiment (Morino et al., 1971): (a) the cross-section z ¼ 0.167 m; (b) the crown line. 6

Y.-X. Peng et al.

Ocean Engineering 192 (2019) 106507

Fig. 10. Deformed configuration with contours of the equivalent plastic strain at load factor P=P0 ¼ 60.

Fig. 13. Load-displacement curves (Inward Load).

1 2 3 62 1 1 0 6 1 6 Q ¼ pffiffi 4 1 1 p0ffiffi 5; ΛP ¼ 6 60 2 0 0 2 4 0

2

3

2

E

3 0

61 ν 0 07 6 7 6 E 7 3 7; ΛCn ¼ 6 6 0 1þν 07 6 6 2 5 4 0 3 0 0

0 0 E 2ð1 þ νÞ

By considering ξnþ1 ¼ QT5 ðσ nþ1 þ βnþ1 Þ, where: � T � 0 Q Q5 ¼ 0 I2 ðkÞ

ðkÞ

7 7 7 7 7 7 7 5

(59)

ðkÞ

(60)

then Eq. (56) can be diagonalized and have a remarkably simple form: (61)

ðkÞ

ξnþ1 ¼ Sξtrial nþ1 Fig. 11. Pinched hemisphere under concentrate force.

where S is a diagonal matrix, and can be expressed as follows: 2

Ξ

1

6 S¼4 0

0

3

7 � ðkÞ � 1 5 γ 1 þ 3Gs nþ1 I2 κ0

(62)

in which:

Fig. 14. Deformed configuration with contours of the equivalent plastic strain at load level F=F0 ¼ 30.

Fig. 12. Load-displacement curves (Outward Load). 7

Y.-X. Peng et al.

Ocean Engineering 192 (2019) 106507

Fig. 15. Design drawing of the ship cabin (Zhu et al., 2004) (a) Cross section; (b) Planform.

�

Table 1 Material parameters of steel.

ðkþ1Þ

1

Material

σy0 /MPa

ν

E/ GPa

Ep /MPa

ρ/(kg m 3)

C/s

Normal steel Low alloy steel

235

0.3

210

250

7800

40.4

5

390

0.3

210

1305

7800

6� 10 4

50.8

Ξ ¼ I3 þ

φ ∂γ φ

ðkÞ

γ nþ1 ¼ γnþ1

� ðkÞ � γnþ1 4 2ΛCn ΛP þ H’ΛP 3 κ0

p

�ðkÞ

(65)

nþ1

The finial stress state is obtained by substituting the convergence solution of γ into Eqs. (52) and (53). In the equation above, we have:

∂γ φ ¼ 2 ξtrial nþ1

�T

⋅ St Ω∂γ S ⋅ ξtrial nþ1

2κ∂γ κ ¼0 κ0

(66)

where ∂γ S ¼ S⋅∂γ S 1 ⋅S, and we have: 2 3 4 0 1 4 2ΛCn ΛP þ H’ΛP 1 5 3 ∂γ S ¼ κ0 0 3Gs I2

(63)

Substituting Eq. (61) into the yield function in k iteration step Eq. (55), we have: ðkÞ � �T T κ2 rnþ1 ðkÞ trial φnþ1 ¼ ξtrial ⋅ S ΩS ⋅ ξ ¼0 (64) nþ1 nþ1 κ0

Considering that κ ¼ κ0

∂γ κ ¼

q and

∂q ∂r

¼

(67) Ep , we have:

∂κ ∂r 2κ ¼ Ep κ0 ∂r ∂γ

The main numerical algorithm is shown in Algorithm 1.

where Ω ¼ QT5 AQ5 is a constant matrix. The trial stress state is invariable during the iterative algorithm, so the yield function φ is a function only related to γ. Then we can solve the variable γ by the Newton iteration:

Algorithm 1.

8

Return mapping algorithm.

(68)

Y.-X. Peng et al.

Ocean Engineering 192 (2019) 106507

Fig. 16. Particle view of the ship cabin while particle spacing choosen as 0.02 m Fig. 18. Particle view of the stiffener; particle spacing: left 20 mm; right 10 mm.

4. Numerical examples

integration cell are adopted Peng et al. (2018), the illustration of the integration cell is shown in Fig. 2. Since the order of the integration scheme is less than the order of the shape function, it shall be regarded as a “reduced integration” by which shear locking is eliminated. Besides, for the integration in the shell thickness, three Gauss points are employed.

To verify the reliability and accuracy of the present meshfree model, we studied several numerical examples in this section. In the RKPM shell formulation, the bilinear shape function and one stress point per

Fig. 17. Contours of displacement of the shipboard at three typical time instants; left: particle spacing is 20 mm; right: particle spacing is 10 mm. 9

Y.-X. Peng et al.

Ocean Engineering 192 (2019) 106507

4.2. Roof loaded by velocity impulse In this section, the present plastic model applied in shell dynamic analysis is considered. A 120∘ cylindrical panel subjected to the velocity impulse is studied, an initial velocity distribution is specified. The problem dimensions and material parameters are summarized in Fig. 5. The curved ends of the roof are simply-supported, and the side bound­ aries are clamped. An initial velocity of 143.51 m/s normal to the shell surface is loaded on a region marked in Fig. 5. Two particle spacing are used to solve this problem. The total number of particles is 556 and 2242 respectively. The displacement histories of the gauging point, initially located at x ¼ 0, y ¼ 0.0762 m, and z ¼ 0.167 m, are shown in Fig. 6. The results calculated by RKPM are compared with the reference simulations using 224 and 4512 Belytschko-Tsay elements respectively (Benson et al., 2010), as well as experimental data from Morino et al. (1971). Since the bilinear shape function with one stress point per integration cell is applied, no locking is observed in the numerical simulation. And it was shown that the re­ sults obtained by the present model are in good agreement with the reference results. Fig. 7(a) gives the residual deformation of the crown line of the roof and Fig. 7(b) gives the residual deformation of a radial cross-section, in which the experimental data is plotted as reference. Fig. 19. Load-displacement curves of the maximum loading point.

4.3. Elastoplastic response of a simply supported plate The inflation of a square plate is presented in this section. This nu­ merical example is popular to verify the shell formulations within finite elasto-plastic strains and has been considered by a number of authors (Eberlein and Wriggers, 1999; Betsch and Stein, 1999; Klinkel et al., 2008). The geometry of the problem and material data are summarized in Fig. 8. The edges of the plate are simply supported that only the vertical displacements u3 are restricted to zero and a constant pressure load P0 ¼ 104 Pa is applied. Here, perfect plasticity without hardening is assumed. In the simulation, the shell is discretized by 49 � 49 particles. The load-displacement curve of the center node of the plate is shown in Fig. 9, and the results obtained by Eberlein and Wriggers (1999) and Betsch and Stein (1999) are adopted as references. In which the load level denotes P=P0 . It was shown that our results are in agreement with the reference results. This example verifies the applicability of the pre­ sent model to the simulation of extreme bending deformations. The deformed configuration with a contour plot of the equivalent plastic strain for a pressure load of P ¼ 60P0 is shown in Fig. 10, which shows a significant change of the shell surface during the deformation.

4.1. A simply supported plate loaded by pressure impulse In this section, the elasto-plastic dynamic response of a simply sup­ ported plate under uniform step pressure is studied. The dimension of the plate and the material parameters are shown in Fig. 2. The initial yield stress is κ0 ¼ 2:068 � 108 Pa, and perfect plasticity without hard­ ening is assumed. In calculation, uniformly distributed particles are used to discretize the structure, the total number of particles are: 5 � 5, 9 � 9 and 27 � 27, respectively. The result calculated by ABAQUS© (64 � 64 S4R elements are used) is adopted as reference. The displacement histories of the center point are shown in Fig. 3, in which the results calculated by the present method is compared with the reference result obtained by ABAQUS©. From which we can see that the results calculated by RKPM converge to the reference results very quickly. The results obtained by 27 � 27 particles are almost identical with the FEM result, which verifies the convergence and accuracy of the present method. Due to the material plastic, the deflection of the center point can not return to its initial state during the vibration, as shown in Fig. 3. In RKPM shell, the average number of nodes in the support domain is about 12, while in the S4R shell element in ABAQUS, the stress point is connected with four nodes. So the calculation efficiency of RKPM is lower than that of FEM. The vibration of the plate in 0.2s is simulated by RKPM and FEM on one Intel i7-7700k processor; the comparison of the computational cost of these two methods is shown in Fig. 4. From which we can see that the time cost by RKPM is about four times than FEM.

4.4. Pinched hemisphere with isotropic hardening A hemisphere loaded by two inward and two outward forces 90� apart. The plastic model with isotropic hardening is adopted. The ge­ ometry parameters and material properties are shown in Fig. 11. This example has been studied previously by a number of authors (Simo and Kennedy, 1992; Eberlein and Wriggers, 1999; Betsch and Stein, 1999). In the simulation, the shell is discretized by 2883 particles. The applied loads F versus the radial displacement is plotted in Fig. 12, and the results presented in previous literatures (Simo and Kennedy, 1992;

Fig. 20. Particle view of the warship. 10

Y.-X. Peng et al.

Ocean Engineering 192 (2019) 106507

Fig. 21. Contours of the Mises equivalent stress of the warship at five typical time instants.

Fig. 22. Depression of the ship hull.

Fig. 23. The deformation of ship at t ¼ 200 ms obtained by RKPM is compared with the result calculated by FEM. 11

Y.-X. Peng et al.

Ocean Engineering 192 (2019) 106507

Fig. 24. Contours of equivalent strain of the warship and the size of the crevasse.

Eberlein and Wriggers, 1999; Betsch and Stein, 1999) are chosen as references. From which, it shows that the results obtained by the present plastic model are in good agreement with the results presented in Eberlein and Wriggers (1999) and Betsch and Stein (1999). Since the stress resultants plasticity model leads to less hardening when compared with J2 model, there is a divergence between the present results and the results obtained by Simo and Kennedy (1992). The deformed configuration of the hemisphere at a load level of F ¼ 30⋅F0 is shown in Fig. 13, in which the contours of the equivalent plastic strains are included (see Fig. 14).

loading. And at the initial moment, due to the reinforcement of the stiffener, the displacement of the outboard shell near the stiffener is relatively small. With the development of deformation, the displacement of the outboard shell reaches the peak. And because of the influence of plasticity, the displacement of the outboard stay at the plateau although the pressure load decreased near to zero. The time history of the displacement at the maximum loading point is shown in Fig. 19, and the result calculated by Ls-Dyna© with mesh size equal to 10 mm is adopted as the reference. From which we can see that the displacement of the gauging point increases rapidly at the beginning and reaches a peak at 1.5 ms. Subsequently, the displacement declines, but due to the plasticity of the material, there is relatively large residual displacement exist and the displacement fluctuates near the final equi­ librium state. By comparing these curves, we can see that the displace­ ment of the gauging point increases while the particle spacing decreases. The calculation converges, because of the displacement of the gauging point almost the same when particle spacing is 10 mm or 8 mm. And the result calculated by RKPM are in good agreement with the FEM result. In Fig. 18, the comparision of the particle views of the stiffener when the particle spacings are 20 mm and 10 mm is given. It can be seen from the figure that there are only three nodes in the width direction of panel and web of the stiffeners when particle spacing is 20 mm, but when particle spacing is 10 mm, there are at least five nodes in the width di­ rection of panel and web of the stiffeners. Considering that the displacement converges when the particle spacing is less than 10 mm, the rule of particle spacing selection can be obtained, that is, at least five nodes are required in the smaller scale direction of the stiffeners.

4.5. Elastoplastic deformation of a ship cabin under impact pressure In the previous sections, the accuracy and effectiveness of present elastoplastic model applied in simple shell structures are verified. In this section, the elastoplastic response of a ship cabin under impact pressure load is studied. The calculation model comes from the research by Zhu et al. (2004), the dimensions of the structure are illustrated in Fig. 15, in which the length unit is millimeter. The defensive longitudinal bulkhead and the main deck are made of low alloy steel and the other components are made of plain steel. The material parameters of these two steel are summarized in Table 1. In the simulation, the linear hardening is adopted and the damage response of the material is ignored. Since the strain rate effect is important in the impact problem, it can not be neglected. In present study, power law expression is adopted. Therefore, the yield stress is obtained: 0 B κðrÞ ¼ @1 þ

�

ε_ peq C

�1p

1 � C A⋅ κ0 þ Ep r

(69)

4.6. Damage response of ship structure under blast shock wave In this section, to verify the applicability of the present meshfree approach in the analysis of real engineering structure, the dynamic response of the warship under the near-field underwater shock wave is studied. In calculation, the average particle spacing is chosen as 0.4 m, in some specialized area the particle spacing may little than 0.4 m, and a total of 123990 particles are used to discretize the structure. The particle view of the warship is illustrated in Fig. 20. In the modeling of ship structure, only the large stiffener and main bulkheads are considered. Besides, the whole stiffened shell structure is discretized by RKPM shell element. Since the pre-processing capability of the present code is still weak, we ignore the small stiffener to simplify the modeling procedure. The thickness of the outboard shell is set to 12 mm while the others shell are set to 7 mm. The material of the calculation model is Low alloy steel, whose material parameters is list in Table 1 in the previous section. In the simulation, 100 kg TNT is placed below the amidship. Considering that the edge length of 100 kg cubic TNT is about 0.4 m, which is equal to the average particle spacing. If we consider the fluidstructure interaction (FSI), the particle spacing should be shrunk at least 4 times, which will result in a tremendous increase in the number of particles. Since the objective is to verify the effectiveness of the struc­ tural algorithm, the FSI effect is neglected and only the underwater shock wave is taken into account. Moreover, in the simulation, no

where ε_ peq denote the equivalent plastic strain rate, C and p is the strain rate hardening coefficients. In calculation, the water-resistant longitudinal bulkhead is clamped, and the impact pressure is loaded on the broadside shell which is rep­ resented in Fig. 15 by a red line. The temporal and spatial distribution of the pressure load are governed by the expression (Zamyshlyaev and Yakovlev, 1973) list below: � �1:51 t 0:2714 p ¼ 44:1 � 106 ⋅e2:84�10 5 (70) jjx xc jj where x is the coordinate of a stress point on the loading surface, and xc denotes the center coordinate which is (0.7,0.25,0.75). In practice, four different particle spacings has been adopted. In these four cases, the total number of particles used to discretize the model are 23274, 41166, 96595 and 145275, respectively. The particle view of the ship cabin while particle spacing is 20 mm is shown in Fig. 16. The fixed time increment is adopted in simulation. The models with different particle spacing are calculated, and the displacement contour of the ship cabin is shown in Fig. 17, while the particle spacing on the left is 20 mm and on the right is 10 mm. From the figure, it can be seen that the displacement of the maximum load point increases rapidly after 12

Y.-X. Peng et al.

Ocean Engineering 192 (2019) 106507

essential boundary conditions are applied, and gravity is not considered. The temporal and spatial distribution of the shock wave are calcu­ lated using the empirical formula (Zamyshlyaev and Yakovlev, 1973). When the distance between the load point and the center of TNT is 2e 6 times to the TNT radius, the magnitude of shock wave is calculated by the formula developed by Ming et al. (2016). In calculation, a fixed time step of Δt ¼ 1:0 � 10 5 s is adopted. The contours of the Mises equiva­ lent stress of the warship at several typical time are shown in Fig. 21. From which, it can be seen that a local depression occurred at the bottom of the amidship as soon as the simulation begins. After that, the defor­ mation of the ship further increases, and finally, the depression of the ship hull is expanded to the range of four longitudinal compartments, which can be observed in Fig. 22. The comparison of the deformed configuration of the ship at t ¼ 200ms calculated by the present method and LS-DYNA is shown in Fig. 23. From which we can see that the results calculated by RKPM are in good agreement with FEM to some extent. In the simulation, the computation took about 304 min on one Intel i7-7700k processor while using RKPM. As a contrast, it took about 3034s while using LS-DYNA, which is about six times less than RKPM. The final contour of the equivalent plastic strain of the ship hull is shown in Fig. 24. From which, we can see that there is a 5.0m � 2.67 m rectangular damage hole is finally produced.

Areias, P.M., Belytschko, T., 2006. Analysis of finite strain anisotropic elastoplastic fracture in thin plates and shells. J. Aerosp. Eng. 19 (4), 259–270. Benson, D., Bazilevs, Y., Hsu, M.-C., Hughes, T., 2010. Isogeometric shell analysis: the reissner–mindlin shell. Comput. Methods Appl. Mech. Eng. 199 (5–8), 276–289. Betsch, P., Stein, E., 1999. Numerical implementation of multiplicative elasto-plasticity into assumed strain elements with application to shells at large strains. Comput. Methods Appl. Mech. Eng. 179 (3–4), 215–245. Caleyron, F., Combescure, A., Faucher, V., Potapov, S., 2012. Dynamic simulation of damage-fracture transition in smoothed particles hydrodynamics shells. Int. J. Numer. Methods Eng. 90 (6), 707–738. Chen, J.-S., Pan, C., Wu, C.-T., Liu, W.K., 1996. Reproducing kernel particle methods for large deformation analysis of non-linear structures. Comput. Methods Appl. Mech. Eng. 139 (1–4), 195–227. Cheng, Z., Thacker, J., Pilkey, W., Hollowell, W., Reagan, S., Sieveka, E., 2001. Experiences in reverse-engineering of a finite element automobile crash model. Finite Elem. Anal. Des. 37 (11), 843–860. Eberlein, R., Wriggers, P., 1999. Finite element concepts for finite elastoplastic strains and isotropic stress response in shells: theoretical and computational analysis. Comput. Methods Appl. Mech. Eng. 171 (3–4), 243–279. Günther, F.C., Liu, W.K., 1998. Implementation of boundary conditions for meshless methods. Comput. Methods Appl. Mech. Eng. 163 (1–4), 205–230. Guo, H., Zhuang, X., Rabczuk, T., 2019. A deep collocation method for the bending analysis of Kirchhoff plate. CMC-COMPUTERS MATERIALS & CONTINUA 59 (2), 433–456. Hughes, T.J., Liu, W.K., 1981. Nonlinear finite element analysis of shells: Part i. threedimensional shells. Comput. Methods Appl. Mech. Eng. 26 (3), 331–362. Karaoulanis, F.E., 2013. Implicit numerical integration of nonsmooth multisurface yield criteria in the principal stress space. Arch. Comput. Methods Eng. 20 (3), 263–308. Kiendl, J., Bletzinger, K.-U., Linhard, J., Wüchner, R., 2009. Isogeometric shell analysis with Kirchhoff–love elements. Comput. Methods Appl. Mech. Eng. 198 (49–52), 3902–3914. Klinkel, S., Gruttmann, F., Wagner, W., 2008. A mixed shell formulation accounting for thickness strains and finite strain 3d material models. Int. J. Numer. Methods Eng. 74 (6), 945–970. Kuipers, J.B., 1999. Quaternions and Rotation Sequences, vol. 66. Princeton university press, Princeton. Li, S., Hao, W., Liu, W.K., 2000. Numerical simulations of large deformation of thin shell structures using meshfree methods. Comput. Mech. 25 (2–3), 102–116. Li, S., Liu, W.K., 2002. Meshfree and particle methods and their applications. Appl. Mech. Rev. 55 (1), 1–34. Liang, C.-C., Tai, Y.-S., 2006. Shock responses of a surface ship subjected to noncontact underwater explosions. Ocean. Eng. 33 (5–6), 748–772. Ming, F., Zhang, A., Xue, Y., Wang, S., 2016. Damage characteristics of ship structures subjected to shockwaves of underwater contact explosions. Ocean. Eng. 117, 359–382. Morino, L., Leech, J., Witmer, E., 1971. An improved numerical calculation technique for large elastic-plastic transient deformations of thin shells: Part 2—evaluation and applications. J. Appl. Mech. 38 (2), 429–436. Nguyen, V.P., Rabczuk, T., Bordas, S., Duflot, M., 2008. Meshless methods: a review and computer implementation aspects. Math. Comput. Simulat. 79 (3), 763–813. Nguyen-Thanh, N., Valizadeh, N., Nguyen, M., Nguyen-Xuan, H., Zhuang, X., Areias, P., Zi, G., Bazilevs, Y., De Lorenzis, L., Rabczuk, T., 2015. An extended isogeometric thin shell analysis based on Kirchhoff–love theory. Comput. Methods Appl. Mech. Eng. 284, 265–291. Nguyen-Thanh, N., Zhou, K., Zhuang, X., Areias, P., Nguyen-Xuan, H., Bazilevs, Y., Rabczuk, T., 2017. Isogeometric analysis of large-deformation thin shells using rhtsplines for multiple-patch coupling. Comput. Methods Appl. Mech. Eng. 316, 1157–1178. O’Grady, J., Foster, J., 2014. Peridynamic plates and flat shells: a non-ordinary, statebased model. Int. J. Solids Struct. 51 (25–26), 4572–4579. Peng, Y., Zhang, A., Li, S., Ming, F., 2019. A beam formulation based on rkpm for the dynamic analysis of stiffened shell structures. Comput. Mech. 63 (1), 35–48. Peng, Y., Zhang, A., Ming, F., 2018. A thick shell model based on reproducing kernel particle method and its application in geometrically nonlinear analysis. Comput. Mech. 62 (3), 309–321. Rabczuk, T., Areias, P., Belytschko, T., 2007. A meshfree thin shell method for non-linear dynamic fracture. Int. J. Numer. Methods Eng. 72 (5), 524–548. Rabczuk, T., Belytschko, T., 2004. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. Int. J. Numer. Methods Eng. 61 (13), 2316–2343. Rabczuk, T., Gracie, R., Song, J.-H., Belytschko, T., 2010. Immersed particle method for fluid–structure interaction. Int. J. Numer. Methods Eng. 81 (1), 48–71. Ren, B., Li, S., 2012. Modeling and simulation of large-scale ductile fracture in plates and shells. Int. J. Solids Struct. 49 (18), 2373–2393. Simo, J.C., Kennedy, J.G., 1992. On a stress resultant geometrically exact shell model. part v. nonlinear plasticity: formulation and integration algorithms. Comput. Methods Appl. Mech. Eng. 96 (2), 133–171. Simonsen, B.C., 1997. Mechanics of Ship Grounding. Ph.D. thesis, Technical University of Denmark (DTU). Truesdell, C., Noll, W., 2004. The non-linear field theories of mechanics. In: The Nonlinear Field Theories of Mechanics. Springer, pp. 1–579. Wang, D., Chen, J.-S., 2004. Locking-free stabilized conforming nodal integration for meshfree mindlin–reissner plate formulation. Comput. Methods Appl. Mech. Eng. 193 (12–14), 1065–1083.

5. Conclusions In present work, a meshfree approach for the numerical simulation of the elasto-plasticity response of the real ship structure is established. The RKPM shell formulation based on Mindlin shell theory is adopted, and a general elasto-plasticity constitutive model for shell structures is pre­ sented. The yield criterion is derived directly through the classic Mises yield criterion, and the zero normal stress assumption is implemented by removing the normal stress component from the yield function. In calculation, the yield function is put into matrix format, by diagonal­ izing the matrices and avoid the stress resultants, the complexity of the return mapping algorithm is reduced. First, several numerical examples are presented to verify the reli­ ability of the proposed elastoplasticity model, the results obtained by RKPM are in good agreement with the reference results, which dem­ onstrates that the present numerical model has good accuracy. Then, the dynamic responese of a ship cabin under impact load is studied, the applicability of the present meshfree approach in dynamic analysis of complex shell structure is verified. And in the simulation, to ensure a reliable result, there should be at least 5 particles in the small-scale di­ rection of the ribs. Finally, the damage response of a real ship structure subjected to underwater explosion shock wave is studied, the capability of the present meshfree approach applied in the engineering problem is verified. Yet, few researches have been done on the application of the meshfree methods in dynamic analysis of complex structure analysis, we believe that the present work lays a solid foundation for future appli­ cation of the meshfree method in complex engineering structures analysis. Acknowledgements The authors thank the National Numerical Wind Tunnel Project of China (2018-ZT2B05) and Industrial Technology Development Program of China (JCKY2018604C010, JCKY2017604C002) for their support. References Amiri, F., Mill� an, D., Shen, Y., Rabczuk, T., Arroyo, M., 2014. Phase-field modeling of fracture in linear thin shells. Theor. Appl. Fract. Mech. 69, 102–109. Areias, P., Rabczuk, T., 2013. Finite strain fracture of plates and shells with configurational forces and edge rotations. Int. J. Numer. Methods Eng. 94 (12), 1099–1122.

13

Y.-X. Peng et al.

Ocean Engineering 192 (2019) 106507

Zamyshlyaev, B.V., Yakovlev, Y.S., 1973. Dynamic loads in underwater explosion. Naval Intelligence Support Center Washington Dc Translation Div, pp. 86–120. Zeng, Q., Combescure, A., Arnaudeau, F., 2001. An efficient plasticity algorithm for shell elements application to metal forming simulation. Comput. Struct. 79 (16), 1525–1540.

Zhu, X., Zhang, Z.-H., Liu, R.-Q., Zhu, Y.-X., 2004. Experimental study on the explosion resistance of cabin near shipboard of surface warship subjected to underwater contact explosion(in Chinese). Explos. Shock Waves 24 (2), 133–139.

14