Numerical Study on Ductile Failure Behaviours of Steel Structures under Quasi-Static Punch Loading

Abstract

A reliable finite element procedure to simulate shear-dominated ductile fractures in large-scale, thin-walled steel structures is still evolving primarily due to the challenges in determining the failure criterion of metal materials under complex stress states. This paper aims to examine the accuracy of the modified Gurson–Tvergaard–Needleman (GTN) model considering the shear failure in simulating the ductile fracture of steel plate structures under quasi-static punch loading. The modified GTN damage models are performed by the ABAQUS user-defined material subroutine (VUMAT). The void-related parameters and shear damage parameters of Xue’s and the N-H modified GTN models are calibrated from test specimens with various geometries corresponding to different stress triaxiality and shearing conditions. The damage evolution associated with shearing of voids in the modified GTN models has strong influences on the stress triaxiality versus plastic strain under complex stress states, especially for the shear-dominated loading conditions. Based on the original GTN model, Xue’s and the N-H modified GTN model with calibrated material parameters, a numerical comparative study examines the ductile fracture of steel non-stiffened plates and stiffened plates under punch loading. Benchmarked against the experimental studies, the numerical simulations demonstrate that the shear-driven void evolution in the modified GTN model imposes significant effects on the load–displacement responses as well as the onset and extension of ductile fractures in steel plates under punch actions. The N-H modified model with calibrated shear damage parameters shows a better correlation with the ductile fractures in steel plates observed in the experiment than the original GTN model and Xue’s modified GTN model. As a result of this study, the modified GTN model considering shear action can be applied for practical applications in the crashworthiness assessment of ship collision and grounding.

1. Introduction

Ship collision and grounding accidents have become common events with increasing tonnage and busier traffic routes. In a collision or grounding accident, the penetration and damage of the ship structure may cause loss of life and property as well as potential pollution [1]. In most cases, ship collision or grounding events are the scenarios of structures subjected to local quasi-static loads with a low loading speed and large loading mass. The crashworthiness analysis of a structure subjected to localised quasi-static punch load is widely used in the safety assessment in the early design stage of a ship. It is very important to develop an effective and accurate assessment procedure for the crashworthiness of a structure under localised load. Consequently, a comprehensive understanding of the failure mechanisms, including the fracture failure driven by shear actions, in marine structures under local quasi-static loads thus becomes essential to enhance the integrity of these structures.

The explicit finite element method is often used to analyse accidental responses of large-scale and thin-walled steel structures due to the rapid development of computational technology. However, the simulation of cracking evolution in marine structures is a challenging task in the ship collision or grounding analysis, especially in the determination of the fracture criterion under large shear actions. Numerical simulations for the penetration responses of a structure including highly nonlinear elastic-plastic deformations until rupture are also intriguing. The elastic-plastic behaviour until the stage of necking can be captured well with the FE method, but beyond the onset of necking, the correspondence to reality for structure fracture may be very poor [2]. Simulations of fracture processes of structures need both accurate material models and failure criteria to address the failure elements. For steel materials under predominantly tensile loading, many researchers have proposed different failure criteria for ductile fracture, including the critical equivalent plastic strain (EPS) criterion, the Rice–Tracey–Cockcroft–Latham (RTCL) criterion [3,4], the GTN damage criterion [5,6,7,8], the Bressan–Williams–Hill (BWH) instability criterion [9], and the Forming Limit Diagram (FLD) [10,11]. Although the BWH criterion and the RTCL damage criterion can predict the resistance versus penetration depth of steel structures under punching action satisfactorily with calibrated material models and appropriate mesh size [12], they cannot be used to reveal the influence mechanisms of material shearing behaviour on the structural responses.

The porous plasticity-inspired micro-mechanical models are widely used in ductile fracture analysis. Some studies [5,8,13,14,15,16,17,18] for ductile fracture analysis have been conducted using the popular GTN model. Gurson [5] first developed a pressure-dependent constitutive relation for yielding a porous material. Tvergaard and Needleman [8] then extended Gurson’s model by improving the yield condition to incorporate the coalescence of micro-voids. Faleskog et al. [13] conducted the micromechanics calibration of void growth considering material strain hardening based on the GTN model. Kim et al. [14] investigated the influences of initial porosity and stress triaxiality on void growth and coalescence. Qian et al. [15,16] investigated the ductile fracture behaviour of tubular joints and ductile crack extension of the API 5L X65 steel pipeline through the GTN plasticity damage model. Gholipour et al. [17] used the GTN model to simulate the ductile fracture of in situ tensile testing specimens under different stress states. Yan et al. [18] carried out numerical simulations for ductile fracture of the cold-formed S700 steel using the GTN model.

Despite its success in simulating ductile fractures in metallic materials under medium to high stress triaxialities, the critical weakness of the GTN model lies in its inability in describing the fracture driven predominantly by shear actions under low stress triaxiality. In pure shear conditions where the stress triaxiality is zero, the void growth behaviour in the GTN model prohibits the void growth, which requires positive mean stress. To deal with this problem, some researchers have modified the GTN model according to the failure criterion of voids in the internal shear band proposed by McClintock et al. [19]. Xue [20,21] developed an extended stress-state function to distinguish tension and shear stress states and modified the GTN model considering the shear actions of voids. Then, Nahshon and Hutchinson [22] proposed a damage law in the GTN model to account for the damage accumulation under shearing. Butcher et al. [23] numerically simulated steel straight tube hydroforming (DP600) using a modified GTN constitutive model to account for the influence of void shape and shear on coalescence. Zhou et al. [24] coupled the volumetric damage and the shear damage into the yield function and flow potential of the GTN model, which can simulate shear-induced damage by combining the damage mechanics method with the void growth model. Malcher et al. [25] assessed two different shear mechanisms proposed by Xue [20] and Nahshon and Hutchinson [22] for the prediction of crack formation. Li et al. [26] developed a modified GTN-Thomason criterion for modelling ductile failures of a metal material. He et al. [27] proposed an ISMGTN model considering shear failure and containing two independent damage mechanisms. This model was used for ductile fracture prediction of aluminium alloys under different stress states.

The literature review shows that the modified GTN models considering the shear influences on coalescence are widely used in ductile fracture analysis under different stress states for some small-scale material specimens. Very little work for ductile failure analysis of a large-scale thin-walled structure using the above-mentioned GTN models is reported. The structures under quasi-static punch load exhibited strong shear failure behaviours in both the accidental scenarios and model tests, as mentioned in Ref. [28]. It is difficult for the current failure models to clearly investigate the influences of material shear failure characteristics on structural responses under quasi-static punch loading, such as the EPS criterion, the RTCL criterion, and the BWH instability criterion. Developing a reliable and accurate numerical failure model to simulate the shear-dominated ductile fracture of large-scale, thin-walled steel structures is still a research difficulty due primarily to the challenges in determining the fracture criterion of steel materials under complex stress states.

Therefore, this paper aims to apply the modified GTN model considering shear effects to numerical failure analysis of steel structures under quasi-static punch loading and assess the influences of the shear damage on coalescence on the resistance to penetration response of the structure. The accuracy of simulations for the ductile fracture of steel structures under quasi-static punch loading using the modified GTN model are verified by experimental results in the paper. This paper starts with the description of the modified GTN model considering shear effects proposed by Xue [20] and Nahshon and Hutchinson [22]. The modified GTN damage models presented in this section are implemented in ABAQUS software via a VUMAT subroutine. The next section presents a series of numerical simulations on the ductile fracture of in situ tensile tests under different stress states using the modified GTN model. In this section, the void-related and damage parameters in the model are calibrated based on the experimental results. The subsequent section carried out numerical simulations for fracture behaviours of non-stiffened and stiffened steel plates under quasi-static punch loading using the modified GTN model. The numerically obtained force–displacement results are compared with the experimental ones. The last section summarizes the conclusions drawn from this study. The research of this paper can provide a useful tool to simulate the failure behaviour of structures in ship collision and grounding.

2. Description of the Constitutive Model

2.1. GTN Model

Gurson [5] initially developed a pressure-dependent yield function from an isolated spherical void in a continuum medium. Then, Tvergaard and Needleman [6,7,8] enhanced Gurson’s model to construct the GTN model theory containing porous plasticity. In the GTN model [5,6,7,8], the yield function of the GTN model is written as

$$\varphi ={\left(\frac{{\sigma }_{eq}}{{\sigma }_m}\right)}^2+2q_1{f}^{*}\cosh \left(\frac{3q_2{\sigma }_h}{2{\sigma }_m}\right)-\left(1+q_3{f}^{*}{}^2\right)=0$$

(1)

where $\varphi$ is the yield function, ${\sigma }_{eq}$ is the von Mises stress of material, ${\sigma }_{eq}=\sqrt{\frac{1}{2}\left[{\left({\sigma }_1-{\sigma }_2\right)}^2+{\left({\sigma }_1-{\sigma }_3\right)}^2+{\left({\sigma }_2-{\sigma }_3\right)}^2\right]}$, and ${\sigma }_m$ is the yield stress as a function of equivalent plastic strain in the matrix ${\overline{\epsilon }}_m^p$; ${\sigma }_h$ is the macroscopic hydrostatic pressure, ${\sigma }_h=p=\frac{1}{3}\left({\sigma }_1+{\sigma }_2+{\sigma }_3\right)$, q₁, q₂, q₃ are the coefficients, and $q_3=q_1^2$. The equivalent void volume fraction f* is a function of the current void volume fraction f:

$$f^{*}=\left\{\begin{array}{l}f\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} f<f_c\\ f_c+\left(\frac{1}{q_1}-f_c\right)\frac{f_c-f_c}{f_f-f_c}\hspace{1em}\hspace{1em}{f}_c\le f<f_f\end{array}\right.$$

(2)

where f_c and f_f are the critical void fraction and the void volume fraction at fracture, respectively. The void coalescence will occur when the void volume fraction reaches a critical value f_c, and the material will fail when the void volume fraction reaches the value of f_f. Equation (2) works very well for low-stress triaxiality cases but slightly worse for high triaxiality cases [29].

The increase in the void volume fraction under loading comes from two parts: the original void growth and the extended void nucleation. The total change rate of the void volume fraction is written as

$$\dot{f}={\dot{f}}_{growth}+{\dot{f}}_{nucleation}$$

(3)

where ${\dot{f}}_{growth}$ is the change rate of voids owing to growth, and ${\dot{f}}_{nucleation}$ is the change rate of voids owing to the extended void nucleation.

The matrix material is assumed to be plastically incompressible. Therefore, the void growth rate can be expressed as

$${\dot{f}}_{growth}=\left(1-f\right){\dot{\epsilon }}^p:\mathbf{I}$$

(4)

where I is the second-order unit tensor, and ${\dot{\epsilon }}^p$ is the macroscopic plastic strain rate tensor. A normal distribution of void nucleation versus the plastic strain is expressed as

$${\dot{f}}_{nucleation}=A{\dot{\overline{\epsilon }}}_m^p$$

(5)

where ${\dot{\overline{\epsilon }}}_m^p$ is the plastic strain rate, and the void nucleation coefficient A controlled by strain is

$$A=\frac{f_N}{s_N\sqrt{2\pi }}\exp \left[-\frac{1}{2}{\left(\frac{{\overline{\epsilon }}_m^p-{\epsilon }_N}{s_N}\right)}^2\right]{\dot{\overline{\epsilon }}}_m^p$$

(6)

where f_N represents the limit value of volume fraction of all particles with potential for micro-void nucleation, ${\epsilon }_N$ is the average equivalent plastic strain, s_N is the standard deviation of the distribution.

In porous materials, the plastic flow is related to the void volume fraction and the accumulated plastic strain. The matrix plastic strain is related to the macroscopic strain by the equivalence of the rate of the plastic work, which can be shown as

$$\mathsf{\sigma }:{\dot{\epsilon }}^p=\left(1-f\right){\overline{\sigma }}_m{\dot{\overline{\epsilon }}}_m^p$$

(7)

where $\mathsf{\sigma }$ is Cauchy’s stress tensor. Therefore, the cumulative plastic strain can be written as

$${\dot{\overline{\epsilon }}}_m^p=\frac{\mathsf{\sigma }:{\dot{\epsilon }}^p}{\left(1-f\right){\overline{\sigma }}_m}$$

(8)

If the hydrostatic pressure effect is ignored, the yield function of the GTN model will be

$$\varphi ={\left(\frac{{\sigma }_{eq}}{{\sigma }_m}\right)}^2-{\left(1-q_1{f}^{*}\right)}^2=0$$

(9)

When f = 0, the yield surface is the original Mises yield surface. The presence of a nonzero f value effectively shrinks the original Mises yield surface, as shown Figure 1. According to Equation (1), the yield surface in the ${\sigma }_{eq}$−${\sigma }_h$ plane under different equivalent void volume fractions can be obtained as shown in Figure 1. The damage accumulation forces a strong dependence of the yield surface and the plastic potential on the hydrostatic stress, in contrast to the original Mises yield surface which is independent of the mean stress.

2.2. Modified GTN model for Shear Failure

The original GTN model cannot capture localization and fracture for low triaxialities in shear-dominated deformations. Under shear actions, the distortion of voids and inter-void linking promotes an effective increase in the material’s internal degradation and contributes to material softening.

2.2.1. Xue’s Modified GTN Model

Xue [20] further extended the original GTN model to incorporate the void shearing damage and assumed that the relative position of the void does not change, as shown in Figure 2. With the shear strain increasing, the distance between the free surface and the boundary of the representative volume element decreases. The cell in the initial undeformed configuration and in the deformed configuration is shown in Figure 2.

In the undeformed configuration, the minimal distance between the free surface and the boundary of the cell is

$$a=L/2-R$$

(10)

The shear strain $\gamma$ in the deformed configuration is

$$\gamma =\tan \alpha$$

(11)

where $\alpha$ is the deformation angle seen in Figure 2b. The minimal distance at a deformed configuration is

$$a^{\prime }=a\cos \alpha =a\sqrt{\frac{1}{1+{\gamma }^2}}a$$

(12)

With the logarithm definition of strain, the artificial strain ${\epsilon }_{art}$ is defined as

$${\epsilon }_{art}=\ln \frac{a}{a^{\prime }}=\ln \sqrt{1+{\gamma }^2}$$

(13)

Additionally, the failure shear strain in the shear band ${\epsilon }_{shearband}$ for a small void volume is defined as

$${\epsilon }_{shearband}=L/2R$$

(14)

The damage with respect to the shearing of the void is expressed as

$$D_{shear}={\epsilon }_{art}/{\epsilon }_{shearband}=2R\left(\ln \sqrt{1+{\gamma }^2}\right)/L$$

(15)

Using Taylor’s expansion for $\gamma =\sqrt{3}{\epsilon }_{eq}$, the artificial strain is reduced to the leading terms, which can be expressed as

$${\epsilon }_{art}=\frac{1}{2}{\gamma }^2$$

(16)

Therefore, Equation (17) can be rewritten as

$$D_{shear}\approx \frac{{\gamma }^2}{\sqrt{\pi /f}}=\frac{3}{\sqrt{\pi }}{f}^{1/2}{\epsilon }_{eq}^2$$

(17)

where $f=\pi R^2/L^2$ is the void volume fraction of the cell. Equation (17) can be only used for a 2D problem.

For the 3D case, a spherical void of radius R is at the centre of the representative cell of length L. A similar relation for damage associated with void shearing can be obtained as

$$D_{shear}=\frac{3}{2}{\left(\frac{6}{\pi }\right)}^{1/3}{f}^{1/3}{\epsilon }_{eq}^2$$

(18)

where $f=4\pi R^3/\left(3L^3\right)$ is the void volume fraction for the 3D cell. The evolution of shear damage can be written using the rate form:

$${\dot{D}}_{shear}=q_4{f}^{q_5}{\epsilon }_{eq}{\dot{\epsilon }}_{eq}$$

(19)

where q₄ = 1.69 and q₅ = 1/2 for a 2D cell structure, and q₄ = 3.72 and q₅ = 1/3 for a 3D cell structure. Xue [20] mentioned that for intermediate to large values of f (f > 0.05), Equation (19) becomes less accurate, and the material constant q₄ should be calibrated based on the experimental results.

The rate of the void volume fraction can be expressed by the sum of the void growth and the void nucleation and the void shearing as

$$\dot{D}={\dot{f}}_{nucleation}+{\dot{f}}_{growth}+{\dot{D}}_{shear}$$

(20)

The Lode angle related to the third invariant of the deviatoric stress tensor is an important parameter in the characterization of the stress rate effect on the ductile failure of the material. The function of the Lode angle proposed by Xue [20] is defined as

$$g_0=1-\left|\overline{\theta }\right|$$

(21)

where g₀ represents the Lode angle function, and $\left|\overline{\theta }\right|$ is the normalized Lode angle which assumes values between $-1\le \left|\overline{\theta }\right|\le 1$. The end conditions of $g_0=0$ for the generalized tension (${\theta }_L=-30°$) and $g_0=1$ for the generalized shear (${\theta }_L=0°$) should be satisfied, where the Lode angle ${\theta }_L$ is

$${\theta }_L={\tan }^{-1}\left[\frac{1}{\sqrt{3}}\left(2\chi -1\right)\right]$$

(22)

$$\theta =\frac{6}{\pi }{\theta }_L$$

(23)

In Equation (22), $\chi =\left(s_2-s_3\right)/\left(s_1-s_3\right)$, and s₁, s₂, and s₃ are the principal values of the deviatoric stress tensor.

2.2.2. Nahshon and Hutchinson’s (N-H) Modified GTN Model

Nahshon and Hutchinson [22] developed a function to quantify the shearing damage effects of voids at low triaxiality, which was directly added to the total change rate of the void volume fraction, although it cannot represent a physical accumulation of the porosity. The extended shear failure term assuming linear relations of the porosity and the deviatoric plastic strain rate is written as

$${\dot{D}}_{shear}=k_{s}f\omega \frac{\mathbf{s}:{\dot{\epsilon }}^p}{q}$$

(24)

The void volume fraction rate is expressed by the sum of the growth of voids and shear effects by

$${\dot{D}}_{shear}=\left(1-f\right){\dot{\epsilon }}^p:\mathbf{I}+k_{s}f\omega \frac{\mathbf{s}:{\dot{\epsilon }}^p}{q}$$

(25)

where k_s is a defined damage parameter used for scaling the magnitude of the damage growth rate in shear loading, and $\omega$ is the function of the stress states, which is defined by

$$\omega =1-{\xi }^2=1-{\cos }^2\left(3\theta \right)=1-{\left(\frac{27J_3}{2q^3}\right)}^2$$

(26)

where $\xi$ is the Lode parameter defined as

$$\xi =\frac{27J_3}{2q^3}$$

(27)

where the q is the equivalent stress, $q={\sigma }_{eq}$, and the third invariant of stress J₃ is

$$J_3=\det \left(\mathbf{s}\right)=\left({\sigma }_1-{\sigma }_m\right)\left({\sigma }_2-{\sigma }_m\right)\left({\sigma }_3-{\sigma }_m\right)={\mathrm{s}}_1{\mathrm{s}}_2{\mathrm{s}}_3$$

(28)

In Equation (14), s₁, s₂, s₃ are the principal values of the deviatoric stress tensor.

A Lode angle dependence function was proposed by Nahshon and Hutchinson [22]:

$$g_0=1-{\xi }^2$$

(29)

where $\xi$ is the normalized third invariant, $-1\le \xi \le 1$.

Figure 3 compares the evolution of the Lode angle functions proposed by Nahshon and Hutchinson [22] and Xue [20].

The shear damage evolutions can be rephrased for arbitrary loading conditions as

$${\dot{D}}_{shear}=\left\{\begin{array}{l}{g}_0{q}_4{f}^{q_5}{\epsilon }_{eq}{\dot{\epsilon }}_{eq} \mathrm{for} \mathrm{Xue}’\mathrm{s} \mathrm{model}\\ g_0{k}_s{\dot{\epsilon }}_{eq} \mathrm{for} \mathrm{N}-\mathrm{H} \mathrm{model}\end{array}\right.$$

(30)

where

$$g_0=\left\{\begin{array}{l}1-\left|\overline{\theta }\right| \mathrm{for} \mathrm{Xue}’\mathrm{s} \mathrm{model}\\ 1-{\xi }^2 \mathrm{for} \mathrm{N}-\mathrm{H} \mathrm{model}\end{array}\right.$$

(31)

The modified GTN model can simulate the fracture mode dominated by the shear mechanism under low stress triaxiality and account for the effective damage accumulation caused by shear mechanisms such as rotation and torsion of micro-voids under the complex stress stage. This model can provide some physical and mechanical explanations for the ductile fracture of metal to a certain extent. However, there are some issues that need to be studied and discussed. The parameters and damage variables in this model need to be identified and calibrated through extensive numerical simulations, and the mesh sensitivity of the model needs to be discussed. This model is mainly used in ductile fracture analysis under different stress states for some small-scale material specimens but is rarely used for failure analysis of a large-scale thin-walled structure. The study for the application of this model in failure analysis of large-scale thin-walled structures needs to be conducted.

3. Quasi-Static Tension Tests Correlations

The finite element method is widely used in failure analysis of different types of structures [30,31,32,33]. The present section presents the numerical simulations for quasi-static tensile tests to study the fracture mechanism of specimens with various geometries under different stress states. The modified GTN damage model is implemented in ABAQUS software through using the VUMAT subroutine. The void-related parameters and shear failure parameters of the modified GTN damage model elaborated in Section 2 are calibrated based on the results of tensile tests. The correlated parameters can be used for modelling the ductile fracture of a complex structure under quasi-static punch loadings in the next section.

3.1. Finite Element Models for Tensile Specimens

In order to calibrate the parameters of the modified GTN model, seven in situ quasi-static tensile tests [32,33] were chosen to validate the numerical results. Figure 4 shows the dimensions of seven tensile specimens with different shapes. The tensile tests for the first four specimens as shown in Figure 4a–d were conducted by Kõrgesaar et al. [32]. The tensile tests for the last two specimens as shown in Figure 4e,f were carried out by Zhang et al. [33]. The force–displacement responses as well as the failure initiation and crack propagation corresponding to each specimen were obtained in the tensile tests [32,33].

The engineering and true stress–strain relations of steel materials can be obtained according to the load–displacement curves from the tensile tests ST-1 and ST-2.

Table 1. Material parameters of specimens.

Material Type	Material 1 [32]	Material 2 [33]
Specimen	ST-1, CH, NT-1, SHS	ST-2, NT-2
Thickness (mm)	3	3.15
Density (kg/m3)	7850	7850
Young’s modulus (GPa)	200	207
Poisson’s ratio	0.3	0.3
$\mathrm{Yield} \mathrm{stress} {\sigma }_y$ (MPa)	275	302.8
Strength coefficient k (MPa)	630	690.2
Strain-hardening index n	0.21	0.2
Plateau strain ${\epsilon }_L$	0.012	0.0189

lists two types of material parameters of steel specimens. For material 1 [32], the true stress–strain relationship follows a Swift model:

$${\sigma }_{eq}=\left\{\begin{array}{l}{\sigma }_y \epsilon \le {\epsilon }_L\\ 1-{\xi }^2 \epsilon >{\epsilon }_L\end{array}\right.$$

(32)

In this equation, ${\sigma }_{eq}$ is the equivalent true stress, ${\sigma }_y$ is the yield stress, and k and n are the parameters of strain hardening. The strain hardening is delayed until the plastic strain reaches the plateau strain ${\epsilon }_L$. The parameter ${\epsilon }_0$ is

$${\epsilon }_0={\left({\sigma }_y/k\right)}^{1/n}-{\epsilon }_L$$

(33)

For material 2 [33], the equivalent true stress–strain relationship is expressed by a pow–law relation:

$${\sigma }_{eq}=\left\{\begin{array}{l}{\sigma }_y \epsilon \le {\epsilon }_L \\ k{\epsilon }^n \epsilon >{\epsilon }_L\end{array}\right.$$

(34)

A series of 3D models were built with the eight-node solid element with reduced integration (C3D8R) in the FE analysis, as shown in Figure 5. The refined and coarse meshes are used within and outside the central region of specimens, respectively, in order to better capture fracture behaviours of specimens and to reduce the simulation times.

3.2. Parameter Calibration of the Modified GTN Model

In the original GTN model, nine material parameters need to be determined and calibrated based on experimental and numerical results, including coefficients of the yield function (q₁, q₂ and $q_3=q_1^2$), void nucleation parameters (f₀, f_N, ${\epsilon }_N$, S_N), and void volume fractions of coalescence and fracture (f_c, f_f). In order to ascertain the range of original GTN model material parameters,

Table 2. GTN model parameters in the literature.

Source	Material	q1	q2	q3	${\mathit{\epsilon }}_{\mathit{N}}$	SN	f0	fN	fc	ff
Tvergaard et al. (1984) [8]	-	1.5	1	2.25	0.3	0.1	0	0.04	0.15	0.25
Skallerud and Zhang (1997) [34]	CMn steel	1.5	1	2.25	0.3	0.1	0.0003	0.006	0.026	0.15
Hambli (2001) [35]	Carbon steel	1.5	1	2.25	0.3	0.1	-	0.04	-	-
Rachik et al. (2002) [36]	Steels DD13, X6Cr17	1.5	1	2.25	0.3	0.1	-	0.04	0.1	0.101
Springmann et al. (2005) [37]	Steel	1.5	1	2.25	0.3	0.1	0.001	0.01	0.01	0.15
Lemiale et al. (2009) [38]	Mild steel XES	1.5	1	2.25	0.2	0.1	-	0.04	0.15	0.25
Marouani et al. (2009) [39]	FeSi (3 wt.%) steel	1.5	1	2.25	0.3	0.1	-	0.04	0.11	0.12
Kossakowshi et al. (2012) [40]	S235JR steel	1.91	0.79	3.65	0.3	0.05	0.0017	0.04	0.06	0.6
Kiran et al. (2014) [41]	ASTM A992 steel	1.5	1	2.25	0.45	0.05	0	0.02	0.03	0.5
Zhao et al. (2016) [42]	DP600 steel	1.5	1	2.25	0.2	0.1	0.0008	0.02	0.028	0.09
Gholipour et al. (2019) [17]	SAE 1010 steel	1.5	1	2.25	0.3	0.1	0.00107	0.00716	0.01	0.15

summarizes the GTN model parameters retrieved from a literature review. Most of the references in Table 2 adopt the same quantities of yield function parameters proposed by Tvergaard et al. [8], q₁ = 1, q₂ = 1.5, q₃ = 2.25, ${\epsilon }_N=0.3$, and S_N = 0.1. For the N-H modified GTN model, the parameter k_s for shear damage needs to be determined. For Xue’s modified GTN model, the parameters q₄ and q₅ need to be calibrated.

Table 3. The value range of modified GTN damage model parameters.

Parameter	q1	q2	f0	fN	fc	ff	q4	q5	ks
Range	1–2	0.8–1.2	0–0.01	0.001–0.04	0.01–0.15	0.1–0.6	1.5–4.5	0.198–0.462	1–4

lists the range of values for material parameters in the modified GTN models. Based on the numerical simulations for the tensile test (ST-2), the sensitivity analysis of the material parameters (q₁, q₂, f₀, f_N, f_c, f_f) of the GTN model is shown in Figure 6.

The comparisons of numerically obtained load–displacement curves under different GTN parameters are given in Figure 6. They show that the final tension displacement increases with the increase in the parameter values (f_c, f_f) or with the decrease in the parameter values (q₁, q₂, f₀, f_N). The parameter of initial void volume fraction f₀ affects the load–displacement response, as it causes a shrink in the yield surface of a material as seen in Figure 1. The other parameters (q₁, q₂, f_N, f_c, f_f) have no effect on the load–displacement curve prior to necking. Figure 6 also shows that the parameter of the volume fraction for void nucleation f_N imposes a strong effect on the load–displacement response of specimen under uniaxial tension. Figure 7 shows the sensitivity analysis of the parameters (q₄ and q₅) for shear damage in Xue’s modified GTN model. As the parameter q₄ decreases or the parameter q₅ increases, the final tension displacement will increase. The larger value of q₄ or the smaller value of q₅ increases the shear damage rate ${\dot{D}}_{shear}$, which facilitates the fracture failure in the specimen.

Based on the sensitivity analysis of the modified GTN model parameters,

Table 4. The calibrated input parameters for Xue’s modified GTN damage model.

Material Type	q1	q2	q3	${\mathit{\epsilon }}_{\mathit{N}}$	SN	f0	fN	fc	ff	q4	q5
Material-1	1.2	0.8	1.44	0.3	0.1	0.001	0.005	0.02	0.35	3.72	0.33
Material-2	1.5	1	2.25	0.3	0.1	0.001	0.005	0.04	0.25	3.72	0.33

lists the calibrated material parameters which lead to the closest agreement between the numerically computed and the experimentally obtained load–displacement curves. To account for the mesh size sensitivity of the modified GTN model, FE simulations are conducted with four types of uniform element sizes in the failure zone. For the specimen ST-1, the mesh lengths are chosen as 0.3 mm, 0.5 mm, 1 mm, and 2 mm, which correspond to four different ratios of mesh size over the plate thickness, 0.1, 0.167, 0.333, and 0.667, respectively. For the specimen ST-2, the mesh lengths are chosen as 0.3 mm, 0.5 mm, 1 mm, and 2 mm, which have the corresponding ratios of mesh size over the plate thickness of 0.0952, 0.159, 0.317, and 0.635.

Figure 8 shows the comparisons of engineering stress–strain curves between standard tensile tests (ST-1, ST-2) and numerical simulations with different mesh sizes based on Xue’s modified GTN model. Before the peak force is reached, the simulated engineering stress–strain curves show a reasonable agreement with the experiments for all mesh sizes. A strong mesh sensitivity in the modified GTN model arises at large deformation levels. Differences between the load–displacement curves become more pronounced when the mesh size increases. The displacement at the onset of fracture with fine meshes is obviously smaller than that with medium or coarse meshes. Figure 9 compares the numerically obtained load–displacement curves based on three material models, the original GTN model, N-H modified GTN model, and Xue’s modified GTN model. In Figure 9, the FE models used have a length of 1 mm. The comparisons show that the predicted final tension displacement using the original GTN model is larger than those using the other two modified GTN models, and the shear damage parameter k_s of the N-H modified GTN model has little influence on the load–displacement response in the standard tension state.

Figure 10 compares the numerically calculated load–displacement curves for NT-1, NT-2, CH, and STS specimens using the original GTN model, N-H modified GTN model, and Xue’s modified GTN model with the experimental results. The numerical load–displacement curves of the NT-1, NT-2, and CH specimens using the modified GTN model compared well with the experimental ones, whereas some slight discrepancies exist for the case of the STS specimen. The force–displacement relations of the original GTN model are highly overestimated, especially for the displacement at the maximal tension before the fracture occurs. The shear damage parameter k_s of the N-H modified GTN model has a significant effect on the prediction of tensile results of the NT-1, NT-2, and STS specimens, except the CH specimen. Figure 11 and Figure 12 present the fracture modes predicted from the FE analysis for the NT-2 and STS specimens. For the NT-2 specimen, the location where the crack initiates and propagates in the numerical simulation of Xue’s modified GTN model agrees with the experimental observation. Due to the differences in the failure criterion, the fracture modes obtained from numerical results of the three GTN models show some divergences. The original GTN model fails to predict the shear failure and leads to an incorrect failure mechanism.

For the tensile shear specimen (STS), the numerical crack initiation and propagation based on the N-H modified GTN model show a good agreement with the experimental observation. In Figure 10d, for Xue’s modified GTN model, the fracture of the STS specimen is captured with a large nominal tension displacement under shear-dominated loading conditions. Compared to the N-H modified GTN model, Xue’s modified GTN model behaves more conservatively, as its peak load is much larger than that of the N-H modified GTN model, and it behaves a much slower post-peak load reduction. For the N-H modified GTN model, the STS specimen fractures with a small nominal tension displacement similar to the experimental data under the same loading condition. Moreover, with the increase in the parameter k_s of shear damage, the material fails quickly with a small tension displacement, as the damage rate ${\dot{D}}_{shear}$ associated with void shearing is positively correlated with the parameter k_s.

Figure 13 illustrates the evolution of the stress triaxiality with respect to the equivalent plastic strain for each tensile specimen. The stress triaxiality is defined as the ratio of the macroscopic hydrostatic pressure and equivalent effective stress of material, ${\sigma }_h/{\sigma }_{eq}$. The evolution of the stress triaxiality is evaluated at the element located in the region (marked by the red dot in Figure 13) where the crack initiates. The different stress triaxiality levels for the same material are observed from the different tensile tests with different geometries. As shown in Figure 13, for specimens ST-1, ST-2, CH, and NT-1, the stress triaxialities at low plastic strain levels obtained from three GTN damage models are very close. Moreover, for the standard tensile specimens (ST-1, ST-2), notched specimens (NT-1, NT-2), and central-hole specimen (CH), the stress triaxialities exceed the value of 1/3 corresponding to the uniaxial tension. In contrast, the initial stress triaxialities in the sheared tensile specimen (STS) remain below the uniaxial tension value of 1/3 at a small deformation level. The predicted stress triaxialities versus equivalent plastic strain for all specimens based on three types of GTN models have some differences as the damage evolves in the element with increasing deformation. The damage evolution associated with the shearing of voids has significant influences on the stress triaxiality versus equivalent plastic strain, especially for the shear-tension-dominated loading conditions.

4. Steel-Plated Structures under Quasi-Static Punch Loading

4.1. Simulations of Non-Stiffened Plate Penetration Test

Zhang et al. [33] presented the experimental study on non-stiffened plates under quasi-static punch loading. Figure 14 illustrates the experimental setup and dimensions of indenter and plate. The thickness of the plate is equal to 3.15 mm. In the test, each of the four edges of the plate specimen had clamp supports between two thick steel plates, which were fixed with the strong base. The plate specimen was loaded at a rate of about 10 mm/min at the centre by hydraulic cylinders. The material used for the non-stiffened plate remained the same as Material-2 in Table 1.

The non-stiffened plate penetration test was analysed using the ABAQUS software with the calibrated modified GTN model given in previous sections. The material parameters of the plate specimen (material-2) are presented in Table 1, and the calibrated input parameters for GTN models are given in Table 4. The eight-node solid element with reduced integration (C3D8R) was used in modelling the plate specimen, as shown in Figure 15. The plate specimen in the contact region was meshed with three different characteristic solid element sizes, 1 mm, 3 mm, and 5 mm. All degrees of freedoms for nodes at the boundary edges of the plate were constrained in the numerical model to simulate the fully clamped boundary condition in the penetration test, without modelling the detailed clamping setup. The indenter was modelled with rigid material because no deformations were observed in the test. The surface-to-surface contact (explicit) defines the interaction between the indenter and plate specimen. The friction coefficient in the contact interaction is defined as 0.2 [33]. The part of the plate specimen in the contact region where the crack initiates and propagates is simulated with fine solid elements.

Figure 16 shows experimental load–displacement curve and numerical load–displacement curves computed from different meshes, GTN models, and damage parameters of k_s. In the test, the contact between the plate specimen and the indenter has three stages as the indenter penetrates through the steel plate: (a) hemisphere–plate contact, (b) cylinder–plate contact, and (c) cone–plate contact, which correspond to different shapes of the three segments of the indenter shown in Figure 14. The whole penetration processes of the non-stiffened plate specimen are categorized into two fracture processes: the crack initiation process and crack propagation process. In the first contact stage (hemisphere–plate contact stage), the experimentally measured load decreases sharply when the crack initiates and propagates with the increase in the penetration depth. In the next contact stage (cylinder–plate contact stage), the measured load remains at a very low level with the increase in penetration depth. In the last contact stage, when the cone part of the indenter punches the ruptured plate, the load increases remarkably, as shown in Figure 16. Figure 16a and

Table 5. The comparisons between experimental and numerical results of the non-stiffened plate test (mesh size = 3 mm).

Results	Test	Simulation Value (Original GTN)	Simulation Value (N-H GTN)	Simulation Value (Xue GTN)
Maximal force (kN)	311.6	365.82	319.52	291.45
Prediction error of maximal force	-	17.4%	2.54%	−6.47%
Displacement corresponding to peak load (mm)	110.25	128.31	108.16	103.69
Prediction error of displacement	-	16.38%	−1.90%	−5.95%

compare the numerical load–displacement curves obtained from the original GTN model, N-H modified GTN model, and Xue’s modified GTN model. The comparisons show that the numerically computed load–displacement response based on the N-H modified GTN model agrees most closely with the experimental data for all three stages of contacts. The modified GTN models considering shear effect show better predictions on the peak load and post-peak response in the experiment than the original GTN model. Figure 16b,c present the numerical load–displacement curves with different mesh sizes obtained from the N-H modified GTN model and Xue’s modified GTN model, respectively. The plate specimen in the contact region is meshed with three different characteristic element sizes, 1 mm, 3 mm, and 5 mm. The mesh sizes of plate model have a marginal effect on the numerical load–displacement relations. For Xue’s modified GTN model, the peak force simulated with coarse mesh size (5 mm) is closest to the experimental value among the three meshes, whereas for the N-H modified GTN model, the fine mesh size (1 mm and 3 mm) has best estimation on the experimental response. Figure 16d elucidates the influences of the damage parameter k_s in the N-H modified GTN model, which shows a negligible effect on the load–displacement responses.

Figure 17 compares the experimentally and numerically obtained fracture shapes of the plate specimen. The N-H modified GTN damage model with k_s = 4 was used in the simulation. In the penetration test, the plate specimen experienced plastic deformations, necking, crack initiation, and propagation, which are also observed in the numerical simulation of the N-H modified GTN damage model. At first, necking appears on the plate with the increase in plastic deformations in the contact region, and the plate subsequently experiences fracture accompanied by a rapid load decrease. Several cracks subsequently emerge in the plate around the necking circle with the punch of the indenter, as shown in Figure 17b. The ruptured specimen separates into several petal pieces as the crack propagates caused by the punching of the conical part of the indenter bar as shown in Figure 17d. Figure 18 shows the fracture shapes of the damaged plate specimen obtained from numerical simulations with different GTN models and mesh sizes. Due to the different failure criteria, the fracture shapes predicted by the three GTN models are quite different. The comparisons in Figure 17 and Figure 18 show that the N-H modified GTN model with a fine mesh size has a most reasonable prediction of crack initiation and propagation in a plate under punch loading.

4.2. Simulations of Laser-Welded Stiffened Plate Penetration Test

Kõrgesaar et al. [32] reported a quasi-static penetration test for a laser-welded stiffened plate performed in the Strength of Material laboratory at Aalto University. Figure 19 shows the experimental setup and fracture shapes of stiffened plate specimen under punch loading. The flat-bar stiffeners with a 30 mm web height and 120 mm spacing were welded on the panel by the laser-welding technique. The in-plane dimension of the stiffened panel was 1200 × 1200 mm, and the thickness of plate and stiffeners equalled 3 mm, as illustrated in Figure 20. The edge of the stiffened plate specimen is fully clamped between two supporting thick plates made from steel S355. A hemispherical indenter with a radius of 120 mm was loaded to punch through the centre of the stiffened plate. The material used for the stiffened plate remained the same as material-1 in Table 1.

The simulation method for a stiffened plate penetration test is the same with that of a non-stiffened plate penetration test. The calibrated modified GTN model is implemented in ABAQUS via a user-defined VUMAT material subroutine. Table 1 and Table 4 list the material parameters of the stiffened plate specimen (material-1) and the calibrated input parameters for GTN models, respectively. The eight-node solid element with reduced integration (C3D8R) was employed in modelling the stiffened plate specimen. The stiffened plate specimen in the contact region was meshed with four types of solid element sizes, 3 mm, 4.5 mm, 6 mm, and 12 mm. All degrees of freedom for nodes on the clamped plate edges remained constrained in the numerical model. The surface-to-surface contact defines the interaction between the indenter and plate specimen. The friction coefficient in the contact interaction is defined as 0.23 [32]. The part of panels and stiffeners in the contact region where the crack initiates and propagates were simulated with fine solid elements, as shown in Figure 21.

Figure 22 illustrates the numerically predicted force–displacement curves with different GTN damage models and mesh sizes and compares them with experiment results. The comparisons of different GTN models in Figure 22a and

Table 6. The comparisons between experimental and numerical results of the stiffened plate test (mesh size = 4.5 mm).

Results	Test	Simulation Value (Original GTN)	Simulation Value (N-H GTN)	Simlation Value (Xue GTN)
Maximal force (kN)	533.39	582.31	526.95	524.38
Prediction error of maximal force	-	9.17%	−1.21%	−1.69%
Displacement corresponding to peak load (mm)	174.18	224.40	188.73	186.97
Prediction error of displacement	-	28.83%	8.35%	7.34%

show that the numerically predicted results based on the modified GTN models considering shear failure are closer to experimental results than those of the original GTN model. Figure 22b,c illustrate the mesh size sensitivity of modified GTN damage models. The peak load and the final displacement decrease slightly with the decreasing mesh size of a stiffened plate model. Xue’s and the N-H modified GTN models with a mesh size of 3 mm have agreed most closely with the experimentally measured load–deformation curves among the four element sizes considered. The mesh sensitivity analysis shows the coarse mesh elements may lead to the over-prediction of load resistance and fail to capture the strain concentrations and hence the crack initiation. Figure 22d shows the effect of the shear damage parameter k_s in the N-H modified GTN model on the force–displacement curve. With the increase in k_s, the peak load and final displacement decline slightly. The numerically predicted load–displacement relation based on the N-H modified GTN model with k_s = 4 and a fine mesh size (3 mm) demonstrates good agreement with the experimental measurement.

Figure 23 presents the numerically predicted post-fracture shape of the stiffened plate with different GTN models and mesh sizes and some test photographs from Kõrgesaar et al. [33]. In the test, cracks firstly initiated on the panel and propagated parallel to the central stiffener, followed by the crack initiation and propagation on the panel in the transverse direction across the stiffener, as shown in Figure 23(a1–a3). Figure 23(b1–b3) compares the numerically computed post-fracture shape of stiffened plate from three types of GTN models. No crack occurs on the stiffened panels based on the original GTN model, whereas some small circular cracks emerge around the central stiffener in the contact region for Xue’s and the N-H modified GTN models. The comparisons in Figure 23(b1–b3) show that the shear-driven failure is a key mechanism in the ductile fracture of plates under quasi-static punch loading. The GTN damage model requires a damage criterion for the shearing of voids to accurately reflect the failure modes in these plates. For Xue’s modified GTN model, the mesh sizes of the stiffened plate model have some influence on the fracture shape, and the finer mesh facilitates the fracture in the plate structures, as shown in Figure 23(c1–c3). Figure 24 shows the cracking path of the stiffened plate model captured from numerical simulations based on modified GTN models. For Xue’s modified GTN model, the crack first initiates and propagates on the panels in the transverse direction perpendicular to the stiffener, followed by crack initiation and propagation in the longitudinal direction along the stiffener, as shown in Figure 24a. For the N-H modified GTN model, the crack first initiates on the panels and propagates parallel to the stiffener, with subsequent crack initiation and propagation in the transverse direction across the stiffener. The stiffener finally ruptures completely as shown in Figure 24b. The fracture path of the N-H modified GTN model resembles what is observed in the experiment but differs from that of Xue’s modified GTN model. Hence, the N-H modified GTN model captures the experimentally observed shear-driven fracture initiation and propagation better than Xue’s modified GTN model.

5. Concluding Remarks

This paper aims to identify an appropriate modified GTN model to simulate the shear-driven fracture failure in steel plates under quasi-static punch loading. Consequently, the original GTN model, Xue’s modified GTN model, and the N-H modified model are investigated and compared to quantify the effect of shear failure on the ductile fracture and the load–displacement responses. The material and damage parameters in the different GTN models are calibrated through the laboratory material tests using specimens with varying levels of shearing and stress triaxiality. The sensitivity of a number of important parameters including the material damage parameters and the element sizes is also examined. The calibrated damage parameters for the various GTM models by the laboratory material tests are subsequently used to investigate the large-scale steel plates under punch loading. Some major concluding remarks of this paper are drawn as follows:

(1)

The sensitivity analysis of void-related parameters (q₁, q₂, f₀, f_N, f_c, f_f) and shear damage parameters (q₄, q₅, k_s) indicate different levels of effects on the force–displacement responses for a standard tensile laboratory material test. The parameters (q₁, q₂, f_N, f_c, f_f) have negligible effects on the force–displacement curves before necking occurs, whereas the parameter of the initial void volume fraction f₀ affects the force–displacement relation before necking because a nonzero f₀ leads to a shrink in the material yield surface. The parameters of q₁, q₂, f₀, f_N, and q₄ have a negative correlation and the parameters of f_c, f_f, and q₅ have a positive correlation with the tension displacement at the onset of fracture. The shear damage parameter k_s of the N-H modified GTN model has a very minor effect on the force–displacement response in uniaxial tension.

(2)

The numerical prediction based on Xue’s or the N-H modified GTN models provides closer estimations on the crack initiation and the load–displacement relationships observed from the tensile tests with various geometries and stress triaxialities and shearing conditions. The damage evolution associated with the shearing of voids has strong influences on the stress triaxiality versus equivalent plastic strain under complex stress states, especially for the shear-driven fracture failure. A strong mesh sensitivity to the failure criterion of the modified GTN model exists in the necking and fracture initiation stages of tests. The differences between the load–displacement curves become more pronounced when the mesh size increases.

(3)

The modified GTN model considering the shear effect with calibrated parameters is applied to simulate steel plate structures under quasi-static punch loading. The experimental results captured from non-stiffened plate and stiffened plate penetration tests are compared well with numerical results based on the N-H modified GTN models. The N-H modified model with k_s = 4 provides much better prediction to the experimental observations compared with the original GTN model and Xue’s modified model. The results demonstrated that the evolution law for the shear failure of voids in the modified GTN model is very important for simulating the ductile fracture of structures under punch loading.

The good agreement between the experimental and numerical results shows that the modified GTN model considering the shear failure can accurately simulate cracking evolution and fracture behaviour in large-scale, thin-walled steel structures. The modified GTN model can be further applied for numerical analysis of failure behaviours of structures in ship collision and grounding events. The study of this paper can provide an effective and accurate numerical calculation procedure for the crashworthiness of a structure under localised load. But for the modified GTN models, it is still a difficult issue to give some quantitative expressions or explanations to calculate the void-related and shear damage parameters, such as f₀, f_N, fc, f_f, and so on, which need in-depth research in the future.