Ductile Fracture Behavior of Mild and High-Tensile Strength Shipbuilding Steels

: A comparison is made of the ductility limits of one mild (normal) and two high-tensile strength shipbuilding steels with an emphasis on stress state and loading path dependency. To describe the ductile fracture behavior of the considered steels accurately, an alternative form of ductile fracture prediction model is presented and calibrated. The present fracture model combines the normalized Cockcroft–Latham and maximum shear stress criterion, and is dependent on both stress triaxiality and Lode angle parameter. The calibrations indicate that, depending on the hardening characteristics of the steels, ductile fracture behavior differs considerably with stress state. It is demonstrated that the adopted fracture model is able to predict the ductile fracture initiation in various test specimens with good accuracy and is ﬂexible in addressing the observed differences in the ductile fracture behavior of the considered steel grades.


Introduction
Ductile fracture resistance is an important structural performance requirement for the steel components of naval structures. As discussed by Matic et al. [1], computational tools are capable of implementing a range of elastic-plastic constitutive relations and ductile fracture initiation models for analyzing the structural problems relevant to ship steels that involve considerable strain beyond moderate plasticity, such as an analysis of the plastic zone ahead of the crack tip and, on a larger scale, the plasticity generated due to forming [2,3] or extreme loads. Therefore, there is a demand for ductile fracture models that can be calibrated and implemented easily in finite element analysis but represent the micro-mechanisms of the failure process as accurately as possible.
Fracture initiation in ductile metals depends strongly on the stress state. Earlier studies on micro-void-based fracture mechanisms highlighted the role of the hydrostatic pressure (or the first invariant of the stress tensor) on void growth [4,5]. Experiments on notched round bars have been conducted to achieve a range of axisymmetric tension stress states. The failure initiation results have been correlated through stress triaxiality (hydrostatic stress normalized with equivalent von Mises stress) and equivalent plastic strain [6,7], and micro-mechanism based models by McClintock [4] and Rice and Tracey [5] have been verified, particularly for high-stress triaxialities. A parallel development is the porous plasticity model by Gurson [8], which was later enhanced by Tvergaard [9], and Tvergaard and Needleman [10]. Porous plasticity models incorporate the void volume fraction into the constitutive behavior of the material, which evolves with void nucleation, growth, and coalescence mechanisms in a phenomenological manner [10].
On the other hand, Johnson and Cook [11] proposed a purely empirical model using a stress triaxiality sensitivity factor similar to the Rice and Tracey model and included the strain rate effect and Appl. Sci. 2020, 10, 7034; doi:10.3390/app10207034 www.mdpi.com/journal/applsci thermal softening terms in a multiplicative format. Cockcroft and Latham [12] proposed an alternative model based on the maximum principal tensile stress to predict fracture initiation in bulk metal forming at a low and negative stress triaxiality. Oh et al. [13] modified the Cockcroft-Latham (CL) model slightly by normalizing the maximum principal tensile stress with the von Mises equivalent stress. CL model incorporates both the hydrostatic pressure (stress triaxiality) sensitivity and the deviatoric stress state sensitivity (so-called Lode angle sensitivity) in an implicit manner. Xue [14], who first interpreted the experimental results of Bao and Wierzbicki [15] correctly, reported the strong dependence of ductile fracture on the deviatoric stress state or the third invariant of the deviatoric stress tensor, i.e., the tendency of one of the principal stresses to dominate over the other two. The test results by Bao and Wierzbicki [15] indicated that the equivalent plastic strain to fracture does not decrease monotonically with increasing stress triaxiality, particularly at a low stress triaxiality range bounded between pure shear and uniaxial tension, and between uniaxial tension and equi-biaxial tension. Xue [14] emphasized that a second measure of the stress state, which is the Lode parameter representing the third invariant of the deviatoric stress tensor, is needed to discriminate between axisymmetric and shear-dominated stress states. Some other test results, such as those by Barsoum and Faleskog [16] and Clausing [17], support Xue's arguments. Consequently, Xue [18] proposed a damage plasticity model for ductile fracture, which is dependent on the hydrostatic pressure and the Lode angle. In addition, Xue [19] developed a Gurson-type model that included the void shearing mechanism because of the Lode angle effect. Based on these recent developments, a dozen empirical or phenomenological uncoupled ductile fracture models that incorporate the Lode angle effect and particularly address fracture due to shear-dominated stress states were proposed recently [20][21][22][23][24][25][26][27]. An important set of these models adopts a simple weighted product format, incorporating terms representing different micro-mechanisms of the ductile fracture initiation process.
A comparative of several popular fracture models was made by Park et al. [28] relating to the prediction accuracy for marine high-tensile strength steels. It was noted that the recent models incorporating Lode angle sensitivity explicitly outperformed the relatively old fracture models such as maximum shear stress criterion (MSS) and CL, yielding a very low error of calibration and high accuracy, especially in the shear-fracture-dominated regime. In addition, the recent models exhibit very similar fracture loci. Therefore, it was noted that the assumptions behind the models are of less importance as compared to the mathematical flexibility that can cover a wide range of stress states. The calibrations and applications of one of these advanced fracture models, namely the Hosford-Coulomb model, were presented by the present authors in a series of recent papers for marine structural steels and low-velocity impact problems [29][30][31][32][33][34][35][36]. On the other hand, despite their ease of calibration and good prediction capability for certain ranges of stress states, MSS and CL criteria lack mathematical flexibility, yielding high calibration errors in general. It was noted by Park et al. [28] that the extension of the CL model by combining it with MSS criterion may provide a flexible alternative format, as demonstrated by Gruben et al. [37]. This paper compares the ductility limits of mild and high-tensile strength shipbuilding steels using an uncoupled ductile fracture model, which is postulated by combining the normalized maximum shear stress and maximum principal tensile stress in a weighted sum format. A hybrid experimental-numerical approach is adopted to identify the hardening law and fracture model parameters for the steel grades considered. The predictive capabilities of the presented fracture model are evaluated for marine structural steels.

Model Formulation
The preliminary information characterizing the stress state and adopted plasticity model is given in Appendix A. Following Gruben et al. [37], a ductile fracture model in association with non-porous plasticity is formulated by combining the normalized maximum shear stress criterion (MSS) and maximum principal tensile stress criterion (modified version of CL criterion by Oh et al. [13]) in a weighted sum format where C 1 ∈ [0, 1] is the weighting factor between the normalized maximum principal tensile stress and the normalized maximum shear stress. C 2 is an exponent that represents the stress state-dependency of the ductility. C 3 is a constant that determines the overall ductility of the material. The operator used for the maximum principal stress, i.e., , is the Macaulay bracket, which yields zero for negative values. Using the transformation of stress states given in Equation (A11) the proposed fracture model can be expressed as follows (2) The fracture model can be expressed as the macroscopic equivalent plastic strain at the instant of first localization under strictly proportional loading paths as follows For general cases, i.e., non-proportional loading paths, the fracture indicator framework is adopted. In this case, the strain to fracture under proportional loading serves as a stress-state-dependent weighting function. The fracture indicator, D ∈ [0, 1], is defined as an integration of the weighted equivalent plastic strain increment along the strain path as follows Note that D is not coupled with the material's elasto-plastic behavior, such that it does not cause any irreversible weakening (softening) with increasing plastic deformation. For finite element analysis of the implentation of the model in the software package Abaqus/Explicit, a user-defined material subroutine (VUMAT) is developed. The user subroutine adopts a standard return-mapping algorithm for plasticity calculations, where equivalent plastic strain increment is determined. Then, it is used in Equation (6) together with the proposed fracture model as the weighing factor to calculate the increment of fracture indicator. An element is assumed to have failed and deleted if all integration points satisfy the condition D = 1.

Parameter Sensitivity
Before discussing the sensitivity of the fracture strain to model parameters, the dependencies of the model to the stress triaxiality and the Lode angle parameter need to be introduced. Figure 1 shows the equivalent plastic strain to fracture obtained using the proposed model with the parameters C 1 = 0.85, C 2 = 2.0, C 3 = 0.2. Geometrically, the fracture locus is represented as a surface in the space of (ε p , η,θ). The fracture locus is asymmetric with respect to the Lode angle parameter. However, its dependence on stress triaxiality is monotonic. 3D fracture locus obtained using the proposed model with the parameters Figure 2 gives an example 2D plot of the fracture locus in the plane of (ε p , η) for several fixed values of the Lode angle parameter. The stress triaxiality dependence is a monotonically decreasing one, which is in agreement with the previous stress-triaxiality-dependent models [4,5,11]. At high stress triaxialities, the effects of the Lode angle vanish and all three curves converge to a single one. Figure 3 shows a 2D plot of the fracture locus in the plane (ε p ,θ) for several fixed values of the stress triaxiality. The asymmetric nature of the Lode angle dependency is more evident in this figure. The Lode angle dependency explains the difference in ductility observed for axisymmetric tension-compression and generalized shear state states.
The parameter C 1 regulates the effects of the maximum shear stress and maximum principal tensile stress. For the limiting case of C 1 = 1, the proposed model reduces to the normalized maximum shear-stress criterion, which is not stress-triaxiality-dependent. The other limiting case, C 1 = 0, reduces the proposed model to the normalized CL criterion (the modified CL criterion proposed by [13]), in which case both the Lode angle-and stress triaxiality-dependencies are retained. Figure 4 presents several plots of the proposed fracture model under plane stress conditions considering several C 1 values. The parameter C 1 has no effect on the strain to fracture for in-plane shear (η = 0). Large C 1 values yield a fracture locus with a noticeable hump at uniaxial tension (η = 1/3) and a plane strain tension (η = 1/ √ 3) "valley", i.e., a trough between uniaxial tension and equi-biaxial tension (η = 2/3). This effect is due mainly to the Lode angle dominance on the stress-state dependency of ductility. In addition, stress trixiality sensitivity decreases with increasing C 1 , as the maximum shear stress criterion dominates. Hence, C 1 plays an important in regulating the stress triaxiality and Lode angle competition. Moreover, a stress-triaxiality-dependent cut-off region value for each C 1 value is noticeable in Figure 4. Ductile fracture will not occur below the cut-off value.     Figure 5 presents several plots of the proposed fracture model under plane stress conditions considering several different C 2 values. The overall effect of C 2 is emphasizing the stress state-dependency by reflecting both the stress triaxiality and Lode angle sensitivities. Figure 6 shows the effects of C 3 . C 3 does not alter the stress-state-dependency of the ductility significantly but shifts the fracture locus up and down. Note that the value of C 3 does not correspond to the strain to fracture for any particular stress state. In the general case, the border of the cut-off region can be derived from Equation (5). The strain to fracture approaches infinity when the following condition is met This equation defines a plane in the space of (ε p , η,θ). Figure 7 shows the projection of the cut-off plane to the plane of (η,θ) for C 1 = 0.6. As mentioned earlier, the parameter C 1 controls the shape of the cut-off plane. For C 1 = 0, the cut-off plane is identical to the biaxial compression part of the plane stress curve, −2/3 ≥ η ≥ −1/3. No cut-off plane exists for the limiting case of C 1 = 1.

Experiments
The experimental data for one normal strength and two high-strength steel grades, which are widely used maritime industry, are used to demonstrate the calibration procedure and validate the model. The target steels are: • A grade (KR) normal-strength (mild) steel plate (6 mm); • AH36 grade (ASTM A131) high-strength steel plate (6 mm); • DH36 grade (LR) high-strength steel plate (7 mm).
Cerik et al. [32] published the test results for grade DH36 steel and Park et al. [28] presented the data for AH36. The tests for A grade are reported in the present study. The base plates with the thicknesses given above are procured from three different Korean steel producers. According to the mill test certificates, Table 1 lists the chemical composition of the steel plates considered. The test specimens are cut from the mid-layer of the plates given above using the CNC wire-cutting method perpendicular to the plate rolling direction. Standard flat dog-bone specimens with a 50-mm-long, 12.5-mm-wide and 2-mm-thick gauge section are used to obtain the flow stress before the onset of diffuse necking. The loading speed is 1 mm/min. In general, marine structural steels display planar-isotropic characteristics, as reported by Park et al. [36]. For simplicity, the mechanical properties of the tested steel grades are assumed to be isotropic. Figure 8 presents the measured true stress versus logarithmic plastic strain curves for the three steel grades up to diffuse necking. Uniaxial tension tests are performed on central hole and notched specimens. These tension specimens are also 2-mm thick. The first type of flat specimen, designated as NT20, has a circular cut-out with a notch radius of 20 mm. The gauge section width is 10 mm. The second type of notched specimen, designated as PST, has a smaller notch radius and but larger gauge width. The central hole specimen (CH) features a 6-mm hole at the center. The shear specimen is adopted from the design proposed by Peirs et al. [38]. Its thickness is 4 mm. Figure 9 shows the scantlings of the test specimens.
The tests are performed with a speed of 0.5 mm/min. The reported displacement correspond to the extensometer measurements at 50-mm gauge length. The tests are conducted three times and the repeatability of the results is confirmed. NT20 CH PST SH Figure 9. Test specimen geometries (dimensions in mm).

Hardening Model
The combined Swift-Voce hardening model, following the general format give by Sung et al. [39], is adopted Here, α ∈ [0, 1] is the weighting factor between the Swift and the Voce hardening laws. The Swift hardening law is given as The Voce law can be written as The hardening law parameters {A, ε 0 , and n}, and Voce law parameters {k 0 , Q, and β} are determined separately using the standard tension test data of flat dog-bone specimens. Only the portion up to the diffuse necking is used. Beyond that point, an iterative procedure is applied to determine the weighting factor α. The test data of NT20 specimen are used for this purpose because the location of fracture initiation in this specimen is precisely defined. Cerik et al. [32] explains this procedure and its rationale in detail. Finite element software Abaqus is used in the analysis of the tests. One eighth of the full gage section of the NT20 specimen is modeled using the symmetry conditions. The model is meshed using eight-node solid elements (C3D8R). The gauge section of the NT20 is meshed with ten elements through half of the thickness. The corresponding element edge length is 0.1 mm. This value is determined after a convergence test with regards to equivalent plastic strain in the gauge section. Figure 10 presents the corresponding hardening curves, and Figure 11 shows the simulation and test results for NT20 specimen. The Swift law overestimates the force levels after the onset of necking, whereas the Voce law yields early localization and softening. It is evident the combined Swift-Voce law provides the closest estimate for the response at large strains.  Figure 11. Comparison of the measured and simulated force-displacement curves of the NT20 specimen.

Loading Paths
The loading paths to fracture initiation, i.e., the history of stress-state parameters (η,θ) at the material point where the fracture initiates, are obtained from the numerical analyses. Analysis of the fracture test specimens (CH, PST, SH) are conducted in a manner similar to NT20 specimen analysis. The same mesh size (0.1 mm) is used in all specimens. Abaqus/Explicit is used with the mass-scaling option to achieve reasonable analysis durations. The total simulation time is 0.01 s, which is long enough to keep the inertial forces negligible. This has been determined by comparing the kinetic energy of the system with the total strain energy. In a quasi-static simulation, as a rule of thumb, the ratio of these two should be lower than 5%. The assumed fracture initiation location is the element that has the maximum equivalent plastic strain at the instant where the experimental force-displacement curve exhibits a sharp decrease. The critical element in CH, PST and NT20 specimens is at the center of the specimens. In the SH specimen, it is on the free surface, slightly away from the center of the gauge. Figures 12 and 13 show the predicted force-displacement curves and test results for the A and AH36 grade test specimens, respectively. It is evident that the force-displacement responses are predicted with good accuracy using the calibrated hardening curves employing combining the Swift-Voce law. Figure 14 shows the loading paths for all three steel grades. The loading paths indicate the evolution of stress-state parameters with increasing equivalent plastic strain. If the stress-state parameters are constant until fracture initiation, the loading path is called proportional. For the CH and SH specimens, the loading paths are relatively proportional, i.e., (η,θ) remain almost constant until fracture initiation. However, PST and NT20 specimens display non-proportional loading paths immediately after necking. Towards fracture initiation, the Lode angle parameter decreases and approaches zero. This observation is in close agreement with the proposed model, where the fracture prefers generalized shear stress states.

Determination of Fracture Model Parameters
To determine the optimal set fracture model parameters χ = {C 1 , C 2 , C 3 }, an optimization problem is set. The loading paths obtained from numerical simulations for each test (η i ε p ,θ i ε p ) and the corresponding fracture indicator are used. All four specimens are considered. The cost function is defined as follows Hence, the optimization problem can be written as The optimization problem is solved with the aid of the Matlab toolbox. In addition, the constraints on C 1 , C 2 and C 3 are constrained as positive values. Table 3 presents the optimal set of the model parameters. Figure 15 presents the resulting 3D fracture loci. The black curves plotted over the 3D fracture loci denote the plane stress condition. Figure 16 shows the plane stress fracture loci.
Several conclusions may be drawn from the calibrated fracture loci. Among the three steel grades, the ductility limits of normal-strength steel, grade A, showed a relatively low Lode angle dependence which could be best described by a fracture locus similar to the Cockcroft-Latham model. On the other hand, the fracture locus of AH36 elucidates the Lode angle effect. Interestingly, AH36 shows higher ductility for certain stress states, such as uniaxial tension and equi-biaxial tension, compared to A grade steel. DH36 showed the lowest ductility among the three steel grades. The fracture locus of DH36 exhibited a relatively high stress-triaxiality-dependence but moderate Lode angle effect. A similar conclusion is drawn by [32]. For negative stress triaxialities and plane strain tension, AH36 and DH36 have similar ductility levels. Figure 17 shows the predicted fracture indicator using the calibrated fracture model at the instant of experimental fracture initiation. Here, the prediction error should be read as deviation from unity. It is apparent that the prediction errors are very low, particularly for the case of the SH specimen. It can concluded that the proposed model could describe the experimental data accurately for all steel grades and test specimens.  Figure 15. Calibrated 3D fracture loci of the steel grades considered.

Conclusions
The experimental results from a series of tests on three different marine structural steel sheets are used to calibrate a proposed ductile fracture initiation model, which is formulated as a weighted sum of normalized maximum shear and principal stresses. A hybrid experimental-numerical approach is adopted. As a by-product of the fracture model calibration process, the flow stresses of the three steel grades at large strains are obtained. The experimental results and simulations are compared in terms of the force-displacement curves and the onset of fracture. Based on the results obtained, the following conclusions can be drawn:

•
Ductile fracture behavior is closely associated with hardening behavior in large strains;

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: The stress state is characterized by the two non-dimensional ratios of the invariants or the principal stresses. In this paper, the stress triaxiality, η, and Lode angle parameter,θ, were used as the stress state parameters. The stress triaxiality (−∞ < η < ∞), which represents the first invariant, is obtained by normalizing hydrostatic stress, σ m , with the von Mises equivalent stress,σ = √ 3J 2 The third invariant can be normalized as follows where θ is the Lode angle (0 ≤ θ ≤ π/3). The normalized Lode angle, which is called the Lode angle parameter (−1 ≤θ ≤ 1), is expressed as follows θ = 0 designates generalized shear (e.g., pure shear or plane strain tension), whereasθ = 1 corresponds to axisymmetric tension andθ = −1 represents axisymmetric compression. An alternative measure of the third invariant is the so-called Lode parameter, L where maximum shear stress, τ max , and normal stress, which acts on the same plane as maximum shear stress, σ n , are defined as follows The stress vector can be defined in the Cartesian coordinate system using the principal stresses (σ 1 , σ 2 , σ 3 ). Alternatively, a cylindrical coordinate system may be used based on (σ, η, θ). The transformation from (σ 1 , σ 2 , σ 3 ) to (σ, η,θ) is obtained as follows In the cylindrical coordinate system, the direction of the stress vector can be defined using stress triaxiality and Lode angle parameter. If these two parameters, which define the stress state, does not vary until fracture, a proportional loading path is experienced.
Under plane stress condition (σ 3 = 0), the stress state can be represented by σ 2 /σ 1 . In that case, the following relation holdsθ