Localized Necking Model for Punching Fracture Simulation in Unstiffened and Stiffened Panels

Abstract

The ductile fracture of thin-shell structures was studied here using a localized necking model. The punching experiments for unstiffened and stiffened panels were compared with numerical predictions using a combined ductile fracture and localized necking model using shell elements. The plasticity and fracture model parameters of JIS G3131 SPHC steel were identified by performing calibration experiments on standard flat bars, notched tension, central hole tension, plane strain tension, and shear specimens. The plasticity beyond the onset of necking was modeled using the Swift hardening law. The damage indicator framework with a combined Hosford–Coulomb fracture model and the domain of shell-to-solid equivalence (DSSE) were adopted to characterize the fracture initiation. The model parameters were calibrated based on the loading paths to fracture initiation, which were extracted from a non-linear finite element (FE) analysis. The presented HC–DSSE model was validated using punch tests and was able to predict fracture initiation with good accuracy.

1. Introduction

Predicting fracture behavior is important for assessing the structural integrity of accident limit states, such as collisions, explosions, and large deformation problems. The application of ductile fracture using nonlinear finite element analysis (FEA) has traditionally been applied to industrial processes, such as sheet metal forming, crashworthiness, and extreme loads. The mixed stress/strain-based fracture model, which is transformed into the space of stress triaxiality, Lode angle, and equivalent plastic strain, is widely used in ductile fracture problems. Therefore, this model was previously demonstrated by experiments and simulations with solid elements sufficiently. Recent studies have focused on accurate fracture prediction with shell elements similar to solid element-based simulations.

In the phenomenological model, fracture initiation occurs when the state variables, represented by stress and strain, reach a critical value. The Rice and Tracey model is based on the void growth theory, which is defined as the ratio of hydrostatic pressure to von Mises equivalent stress [1]. The Cockcroft–Latham model considers the maximum principle stress as more important than the increment of the equivalent plastic strain, and Oh et al. [2] modified this model to include the effect of von Mises equivalent stress. Burak and Choung [3] presented a fracture model in terms of the combination of maximum shear stress and the Cockcroft–Latham model. Bai and Wierzbicki [4] proposed the modified Mohr–Coulomb (MMC) model, which is transformed from a stress-based Mohr–Coulomb criterion to a mixed stress/strain-based fracture model. Morh and Marcadet developed the Hosford–Coulomb (HC) model in terms of the Hosford equivalent stress and normal stress acting on the maximum shear plane [5]. A comparative study of several fracture models was calibrated and validated by Park et al. [6]. It was noted that the old fracture models, such as the maximum shear stress criterion, Johnson–Cook, and Cockcroft–Latham, showed lower predictive capabilities in the low and intermediate stress triaxiality regime, while recent models show a similar fracture locus. The MMC and HC models have been widely used and validated in ductile fracture experiments and simulations with solid elements [3,7,8,9,10,11,12]. Park and Mohr [13] developed a combined HC model with localized necking, which is termed as the domain of shell-to-solid equivalence (DSSE). Alsos et al. [14] proposed the Bressan–Williams–Hill (BWH) criterion to predict fracture ignition with shell elements based on Hill’s localized necking criterion, which is a stress-based FLC. Burak et al. [15] derived an extended strain-based version that considers non-proportional loading from the original BWH. Recently, fracture models with shell elements were applied to the simulation of thin-walled structures [7,13,16,17,18]. The HC-DSSE model consists of fewer material constants than those of other models, which are frequently used because of the low number of experiments required to calibrate the material constants.

We aim to characterize the plasticity and fracture behavior of JIS G3131 SPHC, which is a hot-rolled mild steel plate used for general structures. The plasticity was modeled using the Swift hardening law. The HC–DSSE model, which combines the ductile fracture and localized necking-based fracture models, describes the fracture initiation, and a damage indicator framework is adopted to consider the non-proportional loading paths. The hardening and fracture model parameters were obtained through tensile tests for specimens under various fracture modes. To validate the presented model parameters, structural tests and fracture simulations with shell elements were performed for both unstiffened and stiffened panels. The numerical results were compared with the experimental results, and the prediction capabilities of the proposed fracture model are discussed.

2. Plasticity and Fracture Models

2.1. Stress State Parameters

The stress state for an isotropic solid material can be represented by the non-dimensional parameters, stress triaxiality ($\eta$) and Lode angle parameter $\left(\overline{\theta }\right)$, which are expressed in terms of stress tensor invariants:

(1)$I_1=\mathrm{tr}\left[\sigma \right],$

(2)$J_2=\frac{1}{2}s:s,$

(3)$J_3=\det \left[s\right],$

where $I_1$ is the first invariant of the Cauchy stress tensor $\left(\sigma \right)$, and $J_2$ and $J_3$ are the second and third invariants of the deviatoric stress ($s$), respectively. The stress triaxiality and Lode angle parameters are expressed as follows:

(4)$\eta =\frac{I_1}{3\sqrt{3J_2}},$

(5)$\overline{\theta }=1-\frac{2}{\pi }\arccos \left[\frac{3\sqrt{3}}{2}\frac{\left(J_3\right)}{{\left(J_2\right)}^{3/2}}\right]=1-\frac{2}{\pi }\arccos \left[\frac{27}{2}\frac{J_3}{{\overline{\sigma }}^3}\right];$

$\overline{\theta }$ varies between axisymmetric compression $\left(\overline{\theta }=-1.0\right)$ and axisymmetric tension $\left(\overline{\theta }=1.0\right).$ Based on the modified Haigh–Westergaard coordinates $\left\{\overline{\sigma }, \eta ,\overline{\theta }\right\}$, the principal stress $\left({\sigma }_1\ge {\sigma }_2\ge {\sigma }_3 \right)$ can be reconstructed as follows:

(6)$J_3=\det \left[s\right].$

The stress triaxiality and Lode angle parameters are expressed as follows:

(7)${\mathsf{\sigma }}_1=\overline{\sigma }\left(\eta +\frac{2}{3}\cos \left[\frac{\pi }{6}\left(1-\overline{\theta }\right)\right]\right),$

(8)${\mathsf{\sigma }}_2=\overline{\sigma }\left(\eta +\frac{2}{3}\cos \left[\frac{\pi }{6}\left(3+\overline{\theta }\right)\right]\right),$

(9)${\mathsf{\sigma }}_3=\overline{\sigma }\left(\eta -\frac{2}{3}\cos \left[\frac{\pi }{6}\left(1+\overline{\theta }\right)\right]\right).$

The relationship between stress state parameters $\left\{\eta , \overline{\theta }\right\}$ under plane stress is defined as

(10)$\overline{\theta }=1-\frac{2}{\pi }\arccos \left[\frac{27}{2}\eta \left({\eta }^2-\frac{1}{3}\right)\right].$

2.2. Plasticity

Plastic deformation occurs when the stress due to external loading exceeds the yield stress of the material. The isotropic yield criterion $\left(f\right)$ is defined by the von Mises equivalent stress ($\overline{\sigma }$) and deformation resistance ($k$):

(11)$f\left[\sigma ,k\right]=\overline{\sigma }-k\left[{\overline{\epsilon }}_p\right]=0,$

(12)$\overline{\sigma }=\sqrt{3J_2}.$

The plastic potential theory uses an associated flow rule as follows:

(13)$d{\epsilon }_p=\left(d{\overline{\epsilon }}_p\right)\frac{\partial \overline{\sigma }}{\partial \mathsf{\sigma }},$

where $d{\epsilon }_{\mathrm{p}}$ is the plastic-strain increment vector. As shown in Figure 1, $d{\epsilon }_p$, which occurs when the stress $\mathsf{\sigma }$ reaches the yield surface, is parallel to the tangent plane at the stress point. The equivalent plastic strain increment, $d{\overline{\epsilon }}_p$, is taken as the work conjugate as follows:

(14)$\sigma :d{\epsilon }_p=\overline{\sigma }d{\overline{\epsilon }}_p.$

The deformation resistance, $k$, is a function of the equivalent plastic strain:

(15)$k=f\left[{\overline{\epsilon }}_p\right].$

2.3. Fracture Model

2.3.1. Hosford–Coulomb Fracture Model

Mohr and Marcadet [5] proposed a micromechanically motivated HC model to predict ductile fracture initiation under proportional loading. This model is motivated by the Mohr–Coulomb criterion, which is expressed by combining the Hosford equivalent stress, ${\overline{\sigma }}_{HC}$, and normal stress acting on the plane of maximum shear stress:

(16)${\sigma }_{HC}+c\left({\sigma }_1+{\sigma }_3\right)=b,$

where $b$ is the critical value, and $c$ is the friction coefficient in the Mohr–Coulomb criterion. The Hosford equivalent stress is expressed as follows:

(17)${\sigma }_{HC}={\left\{\frac{1}{2}\left({\left({\sigma }_1-{\sigma }_2\right)}^a+{\left({\sigma }_2-{\sigma }_3\right)}^a+{\left({\sigma }_1-{\sigma }_3\right)}^a\right)\right\}}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$a$}\right.}.$

Equation (17) can be transformed to Haigh–Westergaard coordinates, $\left(\overline{\sigma }, \eta , \overline{\theta }\right)$:

(18)${\sigma }_f\left[\eta ,\overline{\theta }\right]=\frac{b}{{\left\{\frac{1}{2}\left({\left(f_1-f_2\right)}^a+{\left(f_2-f_3\right)}^a+{\left(f_1-f_3\right)}^a\right)\right\}}^{\frac{1}{a}}+c\left(2\eta +f_1+f_3\right)}.$

Finally, the fracture strain for proportional loading, ${\overline{\epsilon }}_f^{pr}$, can be expressed in a mixed stress/strain space $\left(\eta ,\overline{\theta }, {\overline{\epsilon }}_f\right)$ using the inverse of the isotropic power-law type hardening law:

(19)${\overline{\epsilon }}_{HC}^{pr}\left[\eta ,\overline{\theta }\right]=b{\left(1+c\right)}^{1/n_{HC}}g\left[\eta ,\overline{\theta }\right],$

(20)$\mathrm{g}\left[\eta ,\overline{\theta }\right]={\left({\left\{\frac{1}{2}\left({\left(f_1-f_2\right)}^a+{\left(f_2-f_3\right)}^a+{\left(f_1-f_3\right)}^a\right)\right\}}^{\frac{1}{a}}+c\left(2\eta +f_1+f_3\right)\right)}^{\frac{1}{n_{HC}}}.$

The HC model has four parameters: $a, b,c,$ and $n_f$. The transformation coefficient, $n_f$, is set to 0.1, for general steel [5].

Fracture models based on equivalent plastic strain are known to be affected by the loading paths for materials and structures. The linear damage accumulation law was adopted for considering the stress history. For non-proportional loading paths, the HC fracture model is written as follows:

(21)$D_{HC}={{\displaystyle \int }}_0^{{\overline{\epsilon }}_f}\frac{{\overline{\epsilon }}_p}{{\overline{\epsilon }}_{HC}^{pr}\left[\eta ,\overline{\theta }\right]}d{\overline{\epsilon }}_p,$

where ${\overline{\epsilon }}_f$ is the fracture strain. Fracture initiation occurs when the damage indicator D reaches unity.

2.3.2. Localized Necking Model

Shell element-based modeling for numerical analysis of thin-walled shell structures is beneficial for reducing time and simplifying modeling. The fracture model based on the equivalent plastic strain using shell elements is less accurate [17]. Pack and Mohr [13] proposed a domain of DSSE model as an extended HC model for shell element application, which is called HC–DSSE. Figure 2 shows the concept of HC–DSSE in principle strain coordinates $\left\{{\epsilon }_1, {\epsilon }_2\right\}$. The fracture criteria in the HC–DSSE model were defined by the ductile fracture model and necking limit criterion for the fracture criterion. The DSSE model, which is valid in the biaxial tension region $\left(\frac{1}{3}\le \eta \le \frac{2}{3}\right)$ under proportional loading, is similar to the HC model:

(22)${\overline{\epsilon }}_{DSSE}^{pr}=b{\left[{\left\{\frac{1}{2}\left({\left(g_1-g_2\right)}^d+g_1^d+g_2^d\right)\right\}}^{\frac{1}{d}}\right]}^{-\frac{1}{p}},$

where the $g_1$ and $g_2$ are the functions of stress triaxiality:

(23)${\mathrm{g}}_1=\frac{3}{2}\eta +\sqrt{\frac{1}{3}-\frac{3}{4}{\eta }^2},$

(24)${\mathrm{g}}_2=\frac{3}{2}\eta -\sqrt{\frac{1}{3}-\frac{3}{4}{\eta }^2}.$

The recommended value of the coefficient $p$ is 0.01 [13]. $b$ is a parameter in the HC model. Therefore, the parameter in the DSSE model is only $d$. The damage indicator in the DSSE model under non-proportional loading is expressed as

(25)$D_{DSSE}={{\displaystyle \int }}_0^{{\overline{\epsilon }}_f}\frac{{\overline{\epsilon }}_p}{{\overline{\epsilon }}_{DSSE}^{pr}\left[\eta ,\overline{\theta }\right]}d{\overline{\epsilon }}_p.$

Localized necking occurs when $D_{DSSE}$ reaches unity.

3. Fracture Model Formulation

3.1. Material and Experiments

Quasi-static tests were conducted to characterize the plasticity and fracture model parameters of 1.9 mm thick JIS G3131 SPHC hot-rolled thin plates. The chemical composition of the mill test certificate is listed in Table 1. All specimens were fabricated along the rolling direction of the sheet material using the laser cutting technique. Standard uniaxial tension tests were conducted for five types of specimen geometries.

• Standard dog-bone specimen with a 12.5 mm gauge width, denoted as DB [19].

• Central hole specimen with a 6 mm central hole and 20 mm gauge width, denoted as CH.

• A notched tension specimen with a 10 mm gauge width and a notch radius of 20 mm, denoted as NT.

• Plane strain tension specimen with a 16 mm gauge width and a notch radius of 5 mm, denoted as PST.

• Shear specimen, adopted from the design proposed by Park et al. [8], denoted as SH.

All specimen geometries are shown in Figure 3. The tension tests were carried out using a 50 mm extensometer and a servo-mechanical loading frame under a constant stroke velocity of 2.0 mm/min at room temperature. The DB specimen was used to obtain the characteristics of the elasto-plastic behavior before the onset of diffuse necking. The stress state of the CH specimen was close to the uniaxial tension ($\eta =1/3,\overline{\theta }=1.0$). NT and PST specimens show fracture mechanisms similar to plane strain tension ($\eta =1/\sqrt{3},\overline{\theta }=0.0$). The SH specimen was designed to obtain the shear fracture ($\eta =0.0,\overline{\theta }=0.0$) in the gauge part, but the stress state changed to a mixed tension/shear state by increasing the loading step.

3.2. Finite Element Model

Figure 4 shows the FE model of all specimens. The specimens were modeled with eight-node solid elements with reduced integration C3D8R in the Abaqus element library. The half-thickness of the specimens in the gauge section was discretized with ten elements, that is, the element size was 0.095 mm. The FE models were symmetric in the thickness, width, and length, except for SH, that is, one-eighth of the specimen geometry was generated. The SH specimen was symmetric only in the thickness. The force displacement was applied at the top. Material density $\left(\rho \right)$, Young’s modulus $\left(E\right)$, and Poisson’s ratio $\left(\nu \right)$ were 7.85 tonnes/m³, 200 GPa, and 0.3, respectively.

3.3. Plasticity Model Parameters Identification

The Swift hardening law with the Lüders plateau was adopted to model the plastic behavior of JIS G3131 SPHC steel:

(26)$k=f\left[{\overline{\epsilon }}_p\right]=\{\begin{array}{c}{\sigma }_0 \mathrm{if}{\overline{\epsilon }}_p\le {\overline{\epsilon }}_{plat}\\ A{\left({\epsilon }_0+{\overline{\epsilon }}_p\right)}^n \mathrm{if}{\overline{\epsilon }}_p>{\overline{\epsilon }}_{plat}\end{array}.$

The DB and NT specimens were used to identify the Swift law parameters $\left\{A, {\epsilon }_0, n\right\}$ and plateau parameters $\left\{{\sigma }_0, {\overline{\epsilon }}_{plat}\right\}$. The true stress-equivalent plastic strain curve before the onset of diffuse necking was obtained from the tensile test of the DB specimen, as shown in Figure 5.

The flow stress for the post-necking range was extrapolated using Swift’s law, which was validated using the NT specimen. The notched specimen showed a larger fracture strain than that of the DB specimen, making it easier to identify the large plastic strain zone. Figure 6 shows a comparison between the experimental and simulation results for the NT specimen. It seems that extrapolated flow stress based on Swift’s law can be used to predict the material behavior beyond diffuse necking. Table 2 lists the hardening parameters, and Figure 7 shows the extrapolated flow stress with the test results.

3.4. Fracture Model Parameters Identifications

The loading path to fracture is defined as the evaluation of the stress state parameters $\left\{\eta , \overline{\theta }\right\}$ according to the increase in the equivalent plastic strain when fracture initiates. The commercial software package Abaqus/Explicit was used to extract the stress state at the element, which showed the highest equivalent plastic strain at the fracture initiation obtained experimentally. Figure 8 shows a comparison of the force-stroke displacement curves between the test and simulation. The numerical results show a higher accuracy in the prediction of the force levels of the tests. Therefore, the plasticity model worked well. The determination of fracture initiation was defined as a sudden drop in force.

Figure 9 shows the variation in the stress state parameters with increasing equivalent plastic strain. The loading path of the CH specimen was almost constant at fracture, which corresponds to the uniaxial tension stress state $\left\{\eta =1/3, \overline{\theta }=1.0\right\}$. The stress state of the PST specimen initially showed plane strain tension $\left\{\eta =2/3, \overline{\theta }=0.0\right\}$, but the Lode angle parameter changed linearly with the increasing equivalent plastic strain. In the case of the shear specimen, two stress state parameters changed to the mixed shear/tension state from pure shear $\left\{\eta =0.0, \overline{\theta }=0.0\right\}$ because the symmetric plane of the gauge part in the shear specimen rotated with increasing the loading step.

An optimization was conducted to identify the HC model parameters $\left\{a,b,c\right\}$. The objective function $f\left(\chi \right)$ is used to minimize the sum of the square of residues of the numerical error, and the design variables were $\mathsf{\chi }=\left\{a,b,c\right\}$. The optimization algorithm with the linear damage accumulation law is as follows:

(27)$\chi =\arg \min \left\{\wedge \left|\chi \right|\right\}\mathrm{with}f\left(\chi \right)={{\displaystyle \sum }}_{i=1}^3{\left({\left({{\displaystyle \int }}_0^{{\overline{\epsilon }}_f}\frac{d{\overline{\epsilon }}_p}{{\overline{\epsilon }}_{HC}^{pr}\left[\eta ,\overline{\theta }\right]}\right)}_i-1\right)}^2,$

where index $i$ corresponds to the number of calibration tests, which is three for CH, PST, and SH specimens. The MATLAB toolbox with a derivative-free algorithm was used to solve the optimization problem. The calibrated HC fracture model parameters are $a=1.5979,b=1.3599,c=0.0012$.

The DSSE parameter $d$ is identified using the isotropic hardening function $k\left[{\overline{\epsilon }}_p\right]$. The equivalent plastic strain at the onset of localized necking under the plane strain tension condition, ${\overline{\epsilon }}_{DSSE}^{PST}$, is obtained based on the Considère instability criterion:

(28)$k\left[{\overline{\epsilon }}_{DSSE}^{PST}\right]=\frac{2}{\sqrt{3}}\frac{dk\left[{\overline{\epsilon }}_{DSSE}^{PST}\right]}{d{\overline{\epsilon }}_p}.$

Then, parameter $d$ is solved using an implicit equation:

(29)${\left(1+2^{d-1}\right)}^{\frac{1}{d}}=\sqrt{3}{\left(\frac{{\overline{\epsilon }}_{DSSE}^{PST}}{b}\right)}^{-p}.$

For the JIS G3131 SPHC steel, $d$ was 1.7027. Figure 10 presents a 3D HC fracture locus, where the red line corresponds to the plane stress condition. For the HC model, the JIS G3131 SPHC steels tend to be more sensitive to the Lode angle parameter than the stress triaxiality. Figure 11 shows the calibrated HC ductile fracture with the DSSE-localized necking envelope under the plane stress condition. The DSSE model is valid only between uniaxial and equibiaxial tension $\left(1/3\le \eta \le 2/3\right)$.

4. Validation of Fracture Model

4.1. Description of Punch Tests

The unstiffened and stiffened panels were made of the same JIS G3131 SPHC steel used for specimens for calibration and were subjected to punch tests. The unstiffened panel was reported by Park and Choung [17]. The unstiffened panel was an 800.0 mm square, and the stiffener, which had thickness = 1.9 mm and height = 25.0 mm, was welded to the center line of the stiffened panel using single-pass tungsten inert gas (TIG) welding. An indenter was made with an ellipse shape with a major axis of 100.0 mm and minor axis of 50.0 mm. The diagrams of the punch specimens and indenter are presented in Figure 12. The test specimens were supported by a jig and fully edge clamped using 28 M24 bolts and a 5 mm clamping plate, that is, the unclamped effective dimensions were 600 mm square panels. Figure 13 shows the setup of the punching tests and the ruptured plate at the end of the tests. The tests were carried out at a quasi-static speed with an indenter speed of 10 mm/min.

4.2. Numerical Simulation

The numerical simulations for the punching tests were carried out using Abaqus/Explicit with user-defined material subroutines VUMAT for the HC–DSSE model. The effective parts of the punching test specimens were meshed with four-node reduced integration shell elements (S4R in the Abaqus element library), while the indenter was modeled with a rigid element (R3D4 in the Abaqus element library) and located at the specimen center. The general contact option in Abaqus/Explicit was adopted to define the contact between the indenter and specimen, including the friction property, and the friction coefficient was set as 0.23. Abaqus uses five integration points through the thickness of an S4R by default. The first section point is located at the bottom surface of the shell element, and the five-section point refers to the top surface. The damage indicator must be calculated at all the integration points. Subroutine vumatXtrArg with a user-defined material subroutine VUMAT was used to control the state variables in each through-thickness integration point. An element is deleted when one in five thickness integration points meets the condition $D_{DSSE}\ge 1.0$ or $D_{HC}\ge 1.0$.

A uniformly sized element with four different element lengths ($l_e=2.5\mathrm{mm}$, $5.0\mathrm{mm}$, $\mathrm{and}10.0\mathrm{mm}$) were used in the finite element models. The FE model was only generated in the effective parts ($600.0\mathrm{mm}\times 600.0\mathrm{mm}$). Figure 14 shows the FE models of all the specimens and the number of nodes and elements. The boundary condition of fully clamping the edge was used. The mass scaling option was used to decrease the computational cost, and the total simulation time was 0.4 s, which was long enough to neglect the dynamic effect.

4.3. Numerical Results

Figure 15 shows a comparison of the force-stroke displacement curves of punch tests for unstiffened and stiffened panels. The force suddenly dropped at the initiation of fracture. In both cases, the displacement at fracture decreased with the mesh size reduction. Table 3 shows the maximum force,$F_{max}$, displacement at fracture,$d_f$, and simulation error. In both cases, the fracture simulation with a mesh size of approximately 2.5 times the thickness of the plate ($1.9\mathrm{mm}$) most accurately predicted fracture initiation from the experiments.

Figure 16 shows a comparison of the experimental post-fracture shape and numerical simulation $\left(l_e=5.0 \mathrm{mm}\right)$. Table 4 presents the experimental and simulated dimensions of the failure surface for an element length = 5.0 mm. It is confirmed that the presented fracture model with shell elements has an effective predictive capability of not only fracture initiation but also the shape of the fracture surface.

5. Conclusions

The applicability of fracture assessment based on shell elements using a ductile fracture model with localized necking was demonstrated through a structural test for unstiffened and stiffened panels. Tensile tests were carried out to calibrate the plasticity and fracture model parameters of JIS G3131 SPHC steel. The results showed that ductile fracture is more sensitive to the Lode angle parameter than stress triaxiality. For the localized necking envelope, JIS G3131 SPHC showed lower ductility under plane stress tension and proportional loading. It was proved that the presented model showed good accuracy in the prediction of fracture initiation for punching tests. Furthermore, fracture simulations could verify the experimentally obtained post-fracture shapes. Local necking occurs in the local area of the material and structure, and dense elements can simulate the exact necking initiation of the necking in the numerical analysis. Although the HC-DSSE model considers the occurrence and flexible fracture displacement of the necking to be highly approximate, it is difficult to know exactly when the actual fracture occurs after the necking. Therefore, by using dense elements, the precise timing of the necking can be considered fracture, underrating the fracture displacement. This effect can be corrected by making the element size sufficiently large; in this study, structural experiments demonstrated that element sizes ranging from approximately two to five times the thickness of the structure can be used. In the future, it will be necessary to address the element sensitivity part in the fracture model by conducting studies, such as correction of fracture strain, according to the size of the element.

Author Contributions

Conceptualization, S.-J.P.; formal analysis, S.-J.P.; investigation, S.-J.P. and K.K.; methodology, S.-J.P.; software, S.-J.P.; validation, S.-J.P.; visualization, S.-J.P. and K.K.; writing—original draft, S.-J.P.; writing—review and editing, K.K. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by a grant from Tongmyong University Innovated University Research Park (I-URP) funded by Busan Metropolitan City, Republic of Korea.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Geometric representation of the associated flow rule in principle stress space.

Figure 2. Schematic curves of ductile fracture and necking limit.

Figure 3. Specimen geometries and dimensions (unit: mm).

Figure 4. FE models.

Figure 5. True stress-equivalent plastic strain curve of JIS G3131 SPHC.

Figure 6. Comparison of experimental and simulated force-displacement curves for NT.

Figure 7. Flow stress beyond the onset of necking of JIS G3131 SPHC.

Figure 8. The comparisons of force-displacement between experimental and numerical results.

Figure 9. Variation of stress state parameters as function of equivalent plastic strain.

Figure 10. HC fracture locus for JIS G3131 SPHC.

Figure 11. Hosford–Coulomb Fracture locus with localized necking envelope for JIS G3131 SPHC.

Figure 12. Dimensions of the specimens and indenter (unit: mm).

Figure 13. Images of the punching test specimens.

Figure 14. FE models of unstiffened and stiffened panels.

Figure 15. Comparison of numerical results with the test results.

Figure 16. Comparison of experimental and numerical post-fracture shapes.

Table 1

Chemical composition of JIS G3131 SPHC steel.

Material	C	Si	Mn	P	S
JIS G3131 SPHC	0.0509	0.02	0.24	0.014	0.0062

Table 2

Swift hardening law parameters of JIS G3131 SPHC.

$A(MPa)$	${\mathbf{\epsilon }}_{\mathbf{0}}$	$\mathit{n}$	${\mathbf{\sigma }}_{\mathbf{0}}(MPa)$	${\overline{\mathbf{\epsilon }}}_{\mathit{p}\mathit{l}\mathit{a}\mathit{t}}$
624.2	0.0006	0.194	312.72	0.02825

Table 3

Comparison of maximum force and displacement at fracture between experimental and numerical results.

Case			${\mathit{d}}_{\mathit{f}}$ (mm)	Error (%)	${\mathit{F}}_{\mathit{m}\mathit{a}\mathit{x}} (kN)$	Error (%)
Unstiffened panel	Test		72.4	-	82.8	-
	Simulation	$l_e=2.5$	65.9	9.0	81.2	1.93
		$l_e=5.0$	69.1	4.6	82.6	0.24
		$l_e=10.0$	86.1	18.9	103.1	24.5
Stiffened panel	Test		71.9	-	97.6	-
	Simulation	$l_e=2.5$	94.5	31.4	84.1	13.8
		$l_e=5.0$	74.1	3.06	98.2	0.6
		$l_e=10.0$	81.3	13.07	107.6	10.25

Table 4

Measurement of failure surface of experimental and numerical results.

Case		d1 (mm)	d2 (mm)	d3 (mm)	d4 (mm)
Unstiffened panel	Test	7.3	5.1	7.1	8.7
	Simulation	7.32	5.2	6.63	8.04
	Error (%)	0.27	1.96	6.62	7.59
Stiffened panel	Test	12.0	9.9	-	-
	Simulation	11.7	9.73	-	-
	Error	2.5	1.72	-	-