Critical Void Volume Fraction Identification Based on Mesoscopic Damage Model for NVA Shipbuilding Steel

NVA mild steel is a commonly used material in the shipbuilding industry. An accurate model for description of this material’s ductile fracture behaviour in numerical simulation is still a challenging task. In this paper, a new method for predicting the critical void volume fraction fc in the Guson-Tvergaard-Needleman (GTN) model is introduced to describe the ductile fracture behaviour of NVA shipbuilding mild steel during ship collision and grounding scenarios. Most of the previous methods for determination of the parameter fc use a converse method, which determines the values of the parameters through comparisons between experimental results and numerical simulation results but with high uncertainty. A new method is proposed based on the Hill, Bressan, and Williams hypothesis, which reduces the uncertainty to a satisfying extent. To accurately describe the stress-strain relationship of materials before and after necking, a combination of the Voce and Swift models is used to describe the material properties of NVA mild steel. A user-defined material subroutine has been developed to enable the application of the new parameter determination method and its implementation in the finite element software LS-DYNA. It is observed that the model can accurately describe structural damage by comparing the numerical simulation results with those of experiments; thus, the results demonstrate the model’s capacity for structural response prediction in ship collision and grounding scenario simulations


Introduction
Although many efforts have been made to avoid the occurrence of maritime accidents, marine accidents are still inevitable, and collision and grounding have caused nearly half of all previous accidents. An effective way to reduce casualties and economic losses caused by such accidents is to take the ship crashworthiness and survivability after accidents into consideration during the ship design phase. Thus, it is necessary to accurately predict the structural deformation and damage due to collision and grounding accidents. Numerical simulation methods have always been considered as useful tools to obtain the structural deformation and damage with low cost (Calle et al. 2017;Prabowo et al. 2017). However, a widely accepted failure criterion using the finite element method still has not been achieved, due to the difficulty to accurately define failure criteria for the nonlinear ductile properties of NVA steel (Veritas 2007), which is commonly used in shipbuilding.
NVA steel has an intricate ductile failure process in ship collision and grounding accidents. The ductile failure process of metallic materials usually consists of three phases, namely micro-void nucleation, growth, and coalescence. To model the plastic flow and fracture of ductile metal, Gurson (1977) proposed a cell model with micro-voids in a finite matrix and deduced the mesoscopic damage model of porous material. Tvergaard and Needleman (Tvergaard 1982;Tvergaard 1981) modified the original model by introducing three additional fitting parameters. Next, the Guson-Tvergaard-Needleman (GTN) model was proposed, and it uses the void volume fraction to reflect the development of tiny defects in the material. Compared with the traditional plastic model, the yield surface of the Guson-Tvergaard-Needleman model takes into account the effect of the hydrostatic stress and decreases with the increase of the void volume fraction, which reflects the material deterioration characteristic with the damage development in the deformation process. Therefore, the GTN damage model can be considered an ideal mechanical model for describing the ductile fracture process of steel materials in simulating ship collision and grounding accidents.
Failure criteria have always been the focus in research of numerical simulation investigations on ship collisions and groundings (Liu et al. 2015a, b). Several crucial criteria have been proposed in previous research. The combined Rice-Tracey (Rice and Tracey 1969) and Cockcroft-Latham (Cockcroft and Latham 1968) criterion (RTCL) proposed by Törnqvist (2003) has been shown to reasonably predict ductile fractures, and one parameter needs to be determined to calibrate the RTCL criterion. The Johnson and Cook (1983) criterion is believed to be useful in many impact problems with its constants determined by several model tests.
Assuming that the material stress-strain relationship can be expressed by the power law expression, the combined BWH Hill (1952) and Bressan and Williams (1983) instability criterion proposed by Alsos et al. (2008) is able to offer a simplified way to estimate the onset of local necking. Germanischer Lloyd supported the criteria proposed by Lehmann and Yu (1998). The porosity model by Gurson is an advanced damage model, but the main disadvantage of this model is that many parameters cannot be determined directly by simple tensile tests (Törnqvist 2003). Tvergaard and Needleman (Tvergaard and Needleman 1984;Needleman and Tvergaard 1984) did a great deal of pioneering research in meso-damage constitutive modelling and application. Based on the GTN model, many scholars conducted further research. Becker et al. (1989) performed a finite element analysis to explore the effect of void shape on damage development in samples under plane-strain tension. Michel and Suquet (1992) studied the constitutive law of nonlinear porous materials. Thomason (1985) established a two-dimensional model and found that void growth is sensitive to the size and spatial arrangement. Kuna and Sun (1996) used threedimensional cell model analyses to determine that the size and arrangement of voids and interactions have an impact on the yield function. Although the GTN damage model has been continuously developed and used in many fields, the accurate and quick determination of each parameter and establishment of a computationally efficient GTN damage model still require further study. Hogström and Ringsberg (2012) studied the properties of NVA materials in ship collision accidents. In this study, three failure criteria (shear, FLD, and FLSD) are used to initiate damage (DI) at the point of necking, which is followed by a bilinear law for damage evolution (DE) up to the point of fracture, and the results show that the shear criterion with DE can best mimic the experimental results. To reflect the material properties and describe the ductile failure process of NVA mild steel in a more accurate way, a new method is proposed in this paper, based on the Hill, Bressan, and Williams hypothesis and a combination of the Voce and Swift models. Swift and Voce equations are commonly used in flow stress models to describe the stressstrain relationship during material deformation. NVA shipbuilding mild steel has good ductility. It is difficult to simulate the mechanical properties accurately using a single Swift or the Voce equation. To solve this problem, a combination of Voce and Swift equations is used to describe the hardening curve in this study; the two equations are used to describe the stressstrain relationship before and after necking. The GTN model is used to replace the von Mises yield criterion in this study. Furthermore, in previous studies (Mahnken 1999;Bonora 1999;Brunet et al. 2005;Croix et al. 2003), the parameters of the GTN model were usually determined by converse method, which determines the values of the parameters through comparisons between experimental results and numerical simulation results but with high uncertainty. The new method proposed in this paper can reduce the uncertainty of the prediction for critical void volume fraction f c to a certain extent. In this method, the moment of necking is determined by calculating the critical yield stress when the material reaches the necking point, and the void volume fraction f at this moment is determined as f c . In addition, the user-defined subroutine is programmed to realize the application of this method in the finite element software LS-DYNA. The numerical simulation results obtained by the GTN damage model with the new method are compared with the experimental outcomes and numerical results from Hogström and Ringsberg (2012) to demonstrate the effectiveness of the method. This proposed method can contribute to the description of the intricate ductile fracture process of NVA mild steel and to the prediction of structural response in a more accurate way for numerical simulations of ship collision and grounding accidental scenarios.

Description of the Constitutive Model
In this section, all of the theoretical foundations and derivations are introduced. Section 2.1 introduces the development of the GTN damage model and its main theoretical formulas. Section 2.2 describes how to use the stress update algorithm to implement the GTN damage model introduced in Section 2.1 in finite element software LS-DYNA. The selection of a flow stress model to describe the hardening behaviour of the matrix material is explained in Section 2.3, and Section 2.4 presents the proposed new method for determining the parameter f c ; i.e. when the stress reaches the critical value, the value of the continuously updated void volume fraction f mentioned in Section 2.2 is the parameter f c . Gurson (1977) proposed a cell model with micropores in a finite matrix, which can describe the influence of micropores on the mechanical behaviour of the material. To depict the interaction between voids more effectively, Tvergaard (1981Tvergaard ( , 1982 introduced the parameters q 1 , q 2 , and q 3 to make the numerical results more coincident with the experimental results. Besides, Tvergaard and Needleman Needleman and Tvergaard 1984) introduced a damage function (f * ) to explain the phenomenon of sharp polymerization of voids during material rupture. Consequently, the GTN model can be expressed as:

GTN Damage Model
where Φ is the yield function, σ eq is the macroscopic von Mises equivalent stress, σ m is the yield stress of the matrix material, and σ h = −σ kk /3 is the hydrostatic stress. Three parameters q 1 , q 2 , and q 3 modify the original model to better represent void interaction effects. Normally, q 1 = 1.5, q 2 = 1.0, and q 3 = q 1 2 . When taking q 1 = q 2 = q 3 = 1.0, the GTN model degenerates to the original Gurson model. f * is the total effective void volume fraction to account for the gradual decrease of material load-carrying capacity due to the void polymerization. f * approaching zero when the material is undamaged, and the yield function φ degenerates into the standard von Mises form. It is defined as: where f c is the critical void volume fraction at which void Þ is the acceleration factor of the void growth, and f f is the void volume fraction at the moment of material fractures. In addition to the void volume fraction f, plastic flow of the porous material also depends on the equivalent plastic strain of the matrix material. Based on the principle of equivalent plastic work, the evolution function can be obtained as: where dε pl m is the equivalent plastic strain increment of the matrix material, and dε p is the plastic strain increment.
The total damage evolution of the material consists of two parts: the existing void growth and new void nucleation. The increment of the void volume fraction can be expressed as: Assuming that the matrix material is incompressible, the existing void growth df growth is related to the hydrostatic component of plastic strain: where I is the second-order unit tensor.
Assuming that the void nucleation is controlled by plastic strain, the change of the void volume fraction caused by nucleation can be expressed as (Chu and Needleman 1980): where A is the void nucleation coefficient controlled by plastic strain: where f N is the void volume fraction of the nucleating particles, ε N is the mean equivalent plastic strain at void nucleation, and s N is the standard deviation of the nucleation strain.

Numerical Implementation of the Constitutive Model
In this section, the numerical integration algorithm of the constitutive equation, i.e. the stress update algorithm, is described in detail. The implementation details of GTN damage model in finite element software LS-DYNA and a flow chart are given in Fig. 1.
As a satisfactory mathematical model describing the ductile fracture process, the GTN model plays an important role in numerical simulation. Combining the backward Euler complete implicit algorithm proposed by Aravas (2010) with the finite element explicit solution, the implicit calculation of large stiffness matrix by finite element solver can be avoided when the stress update of GTN model is realized, which is especially applicable for sheet metal deformation analysis.
The implementation of stress update algorithm usually includes two parts: elastic prediction and plastic correction. The total strain is separated into elastic and plastic parts: ε = ε e + ε p . The plastic strain increment also consists of two parts: ∂σ eq , n is the flow direction, and λ is the scale factor. The void volume fraction f and the matrix equivalent plastic strain ε pl m are taken as two internal variables of the model. f t and f t + Δt represent the value f at time t and t + Δt respectively. The concrete realization steps are summarized as follows: Step 1. Get the initial condition at time moment t t σ; t ε; tþΔt Δε; t H a where t H a represents the void volume fraction and matrix equivalent plastic strain.
Step 2. Calculate the trail elastic stress by assuming that the strain increment is purely elastic where C is the fourth rank elastic modulus. The hydrostatic stress of trial elastic stress is calculated by: tþΔt σ e kk , and the equivalent stress is calculated by Step 3. Calculate the yield function and determine whether plastic deformation occurs If t + Δt Φ e ≤ 0, the material is still in an elastic state at the current time step, and the true stress is the trail elastic stress. Go to step 5 directly. If t + Δt Φ e > 0, the material is in a plastic state, and the plastic correction in step 4 is needed.
Step 4. Calculate the plastic correction To simplify the presentation, neglect the left superscript t + Δt in this step.
Calculate the flow direction according to Eq. (10): The Newton-Raphson method is used to solve the nonlinear equations. For three-dimensional solid elements, the flow rule and yield condition need to be satisfied simultaneously.
where Δε pl h and Δε pl eq denote the increment of hydrostatic plastic strain and equivalent plastic strain, respectively. Iteration continues until |f 1 | and |f 2 | are both less than the tolerances, which signify the iteration converges; otherwise, terminate the calculation.
For the shell elements, in addition to the above conditions, the stress in the thickness direction of the elements should be zero, that is, σ 33 = 0. After rewriting, the additional condition can be expressed as: where Δε 3 is the normal component of the strain increment in the x 3 direction. Step 5. Update stress and statement variables For solid elements, the hydrostatic stress and equivalent stress are: For shell elements, they can be expressed as: ( where K and G represent the bulk modulus and the shear modules, respectively. Update the stress and statement variables: If Step 6. Go to the next time step Figure 1 shows the flow chart diagram for stress update of GTN damage model.

Matrix Hardening Behaviour
The Swift (1952) and Voce (1948) equations are commonly used expressions to describe the material stress-strain curve.
To determine an appropriate hardening model, the uniaxial tensile test performed by Hogström et al. (2009) is used. Specific experimental details and the determination process can refer to Section 3.2.1, here no longer say.
A single Swift equation is verified firstly. The Swift model can be expressed as where σ m is the yield stress,ε pl m denotes the equivalent plastic strain, ε 0 is the yield strain, K is the hardening coefficient, and n is the hardening exponent. The results are as follows: K = 783 MPa, ε 0 = 0.02, n = 0.26, and the comparison of the experimental data and numerical results is shown in Fig. 2. Figure 2 shows that if the von Mises yield criterion or the GTN damage model is adopted, the Swift hardening model can describe the stress-strain relationship of the material accurately during the initial stage of deformation, which shows that the void volume fraction of NVA mild steel increases slowly at the initial stage of deformation before the neck shrinkage occurs. However, no matter which yield criterion is adopted, the single Swift hardening model cannot accurately describe the mechanical behaviour of the NVA shipbuilding material after constriction.
As for the Voce equation, similar conclusions are observed. Therefore, the Voce and Swift formulas are used to describe the stress-strain relationship of NVA shipbuilding steel before and after constriction.

Determination of the Parameter f c
In this section, a new critical void volume fraction identification method is proposed. It has been demonstrated that using a single Voce or Swift equation cannot simulate the mechanical behaviour of NVA mild steel accurately, so the combination of two equations is used to describe the hardening curve of the material in this paper. The stress-strain relationship can be expressed as: where σ y represents the initial yield stress, ε c is the necking strain, and Q, b, K, B, and n are the hardening parameters of Fig. 2 The comparison of the experimental and numerical results calculated by the Swift equation each model. To ensure the continuity of yield stress before and after necking, B can be expressed as: Experiments have shown that the void volume fraction of carbon steel increases slowly before necking occurs (Törnqvist 2003;Roy et al. 1981), so the necking moment of NVA mild steel can be approximately determined based on the von Mises yield criterion. Hill proposed a criterion for local necking in the negative β regime (Hill 1952). He assumed that a local neck will form with an angle tan À1 1= ffiffiffiffiffiffi −β p À Á to the direction of the major principal stress, which yields reasonable results only for the negative values of β. At the moment when a neck is formed, the strain hardening effect and the thickness reduction balance each other. It means that the traction force reaches a maximum value at the moment of local necking, so the increment of the traction is equal to zero, which leads to the following local necking criterion (Alsos et al. 2008) where β represents the strain rate ratio,β ¼ dε 2 dε 1 , dε 1 and dε 2 are the increment of principal strain, and σ 1 denotes the principal stress. Besides, Hill assumed that the ratio of principal stress increments is consistent with the principal stress ratio, i.e. Table 1 Material parameters Considering the situation of plane stress state, the following expressions can be obtained by von Mises yield function where dσ m and dσ represent the increment of yield stress and equivalent stress, respectively. Combined with Eq. (19), Eq. (23) can be rewritten as According to the content of the plastic mechanics, Eq. (24) can be further written as where the equality holds up if and only if the necking occurs. Combined with the stress-strain relationship before the neck is formed, i.e.σ m ¼ σ y þ Q 1−e −bε pl m , the increment ratio of yield stress and equivalent plastic strain can be also expressed by By contrasting Eq. (25) and Eq. (26), the critical yield stress can be obtained when the necking occurs.
Combined with Eq. (22), the critical yield stress can be expressed by As mentioned above, Hill's assumption yields reasonable results only for the negative value of β. In the positive regime, the method of estimating the onset of local necking proposed by Bressan and Williams (Bressan and Williams 1983) is used. The BW criterion has a simple expression and has been widely used, and it can be expressed by where τ cr is the critical shear stress, which can be calibrated either from uniaxial tensile tests or biaxial tests. Another alternative way used in this paper is calibration at plane strain, β = 0, from the analysis based on Hill's assumption. Combined with Eq. (28) at plane strain, the critical shear stress can be expressed by  As the strain rate ratio becomes negative, the effectiveness of BW criterion becomes questionable (Alsos et al. 2008). To cover the full range of β, Eqs. (29) and (30) are combined to determine the moment when a neck is formed With continuous structural deformation, the yield stress increases gradually. Necking occurs when the value of critical yield stress in Eq. (32) is reached. The necking strain ε c can be derived by substituting critical yield stress into Eq. (20), and the critical void volume fraction f c is the void volume fraction f at this moment, which is continuously updated in Eq. (16) and Eq. (17) in Section 2.2. The assumptions and method proposed in this chapter are verified by uniaxial tension tests, and a ship-like structure experiment introduced in the next section; the results demonstrate the feasibility of the proposed method.

Subroutine Checking
Plane stress tension was used to verify the implementation of the original GTN damage model, as illustrated in Fig. 3. The single element type is a shell, and the initial element dimension is taken as 1 mm. The relationship between the yield stress and equivalent plastic strain of the matrix material is described by the Voce equation, . The material parameters selected are listed in Table 1. The uniaxial tensile simulation results of the plane stress element are shown in Fig. 4. The tensile stress versus strain and void volume fraction versus strain obtained by the origin GTN model ('User defined' line in Fig. 4) are compared with those calculated by the 120-GURSON model ('Explicit' line in Fig. 4) in LS-DYNA constitutive model libraries. It shows that the results obtained by two different ways are nearly the same, which validates the correctness of the original GTN model.

Determination of NVA Mild Steel Material Parameters
To describe the mechanical properties and damage behaviour of the material, the uniform deformation stage in the uniaxial tensile test should be used to evaluate the mechanics property parameters, and the GTN damage parameters can be obtained by the comparison of simulation results with those from the experiments. In this section, the parameters of the GTN damage model are determined by a uniaxial tensile test, and the correctness of the identified parameters is further verified in Section 3.3.
(a) Before fracture (b) After fracture Fig. 9 Distribution of void volume fraction f before and after fracture occurs Fig. 10 The comparison of engineering stress-strain curve from numerical simulation and uniaxial tensile test

Determination of Mechanical Property Parameters
Uniaxial tensile tests are commonly used to evaluate the mechanical properties of materials. In this paper, the mechanical property parameters of NVA mild steel are determined according to the tensile test carried out by Hogström et al. (2009). The specimens used in this test were manufactured according to the DNVGL rules, and the geometric model is shown in Fig. 5. The transition radius R of the test specimen is 25 mm, and the plate thickness a and width b are 4 mm and 25 mm, respectively. The gauge length L o is 56.6 mm, and the parallel test length L c is 76.5 mm. Figure 6 shows the engineering stress-strain curve for NVA mild steel from the test. When the maximum value of engineering stress is reached; i.e. the instant a neck is formed, the material deforms uniformly, and then the engineering stress begins to drop until the fracture occurs. As shown in Fig. 6, the NVA mild steel has good ductility.
To obtain the mechanics property parameters of the material, the engineering stress-strain curve needs to be transformed into the true stress-plastic strain curve, in which the uniform deformation phase of the material is fitted with the hardening model. Voce equation is used to describe the hardening curve before necking occurs.
According to the experimental data shown in Fig. 6, the NVA material property parameters can be obtained by using data analysis software OriginPro 9.1. The initial yield stress and hardening parameters are as follows: σ y = 287.7 MPa, Q = 244.9 MPa, and b = 11.10. The fitting result is shown in Fig. 7. It can be seen that a good agreement is obtained between the Voce hardening curve with the selected parameters and the stress-strain curve of NVA mild steel in the uniform deformation stage, and the correlation coefficient reaches 0.995.

Determination of GTN Damage Parameters
Most damage parameters in the GTN model cannot be directly determined by corresponding experiments; therefore, these parameters should be regarded as tuneable parameters. By adjusting these damage parameters to make the numerical simulation results of the uniaxial tension test consistent with the experimental outcomes, the value of the selected parameters can be considered to reflect the damage evolution law of the material. Based on the user-defined subroutine verified in Section 3.1, the new method to determine the void volume fraction f c is implemented in this section. The GTN damage Fig. 12 The engineering stress-strain curves from numerical simulations and test  The finite element model of a tensile specimen is shown in Fig. 8. The model is established by shell elements with four nodes, the number of shell thickness integration points is set to two, and the mesh size is set to 1.9 mm. The stress-strain relationship is expressed by Eq. (20), and the material parameters are shown in Table 2. The strain rate effect is neglected due to the experimental low-speed loading of 0.1 mm/s. One end of the model is rigidly fixed, and displacement loading is applied on the other end. Figure 9 shows the distribution of void volume fraction before and after a fracture occurs. It can be seen that the damage degree of the tensile specimen increases gradually with the deformation development. When the neck is formed, the damage mainly concentrates in the symmetric centre of the sheet metal, and failure occurs at the instant f f is reached.
The engineering stress-strain curve for NVA mild steel from numerical simulation of uniaxial tension is shown in Fig. 10. A good agreement is obtained between the numerical results and experimental data, reflecting the accuracy of selected material parameters and the feasibility of the proposed method of identifying critical void volume fraction.
To explore the relationship between the mesh size and the void volume fraction at the final fracture f f , the finite element models with mesh sizes of 2.9, 4.8, 6.4, and 9.6 mm are built.
The parameter f f is determined inversely by comparing the numerical results and test data, and the f f versus mesh size curve is drawn in Fig. 11. The numerical results of each finite element model compared with the uniaxial tensile test outcomes are illustrated in Fig. 12. It can be observed that the experimental data are well-fitted, which reflects the reasonability of the selected parameters f f .

Numerical Validations and Discussion
To verify the proposed GTN damage model with the new identifying critical void volume fraction method, a numerical simulation for an experiment of the side-shell structure subjected to collision resistance was conducted. Karlsson et al. (2009) conducted the experiment, and the experimental setup is shown in Fig. 13.
A reinforcing frame was designed and welded around the structure along its edges to create clamped boundary conditions to ensure well-controlled failure modes of the structure (Karlsson et al. 2009). The lower part of the frame was welded to a rigid fixture, see Fig. 13. The finite element model is illustrated in Fig. 14, which is established by shell elements with a mesh size of 10 mm. Substantially, all the elements are four-node Belytschko-Tsay elements (BLT), while few elements located at the intersection of components are threenode BLT elements for the continuity of the structures. BLT element type is selected due to its fast calculating speed and its numerical stability for large deformation problems. Considering the computational accuracy and efficiency, the number of shell thickness integration points is set to five. The value of mesh resolution affects the selection of parameter f f and also the structural deformation mode. To capture the buckling collapse and fracture deformation mode precisely with the reasonable computation time, the mesh size is set to be 10 mm for most of the structural members by research. It can be seen that the selected mesh resolution is accurate enough to obtain a satisfying simulation result, and the computation time is also acceptable.
The material of the side-shell structure is defined as NVA mild steel, where the material parameters used are listed in Table 2. It is worth noting that the strain rate effect is disregarded due to the low-speed loading of the experiment.
According to the f f versus mesh size curve drawn in Fig. 11, (a) Before fracture (b) After fracture Fig. 15 The distribution of void volume fraction on the upper plate before and after fracture the parameter f f of this model is set to 0.05; i.e. the finite element will fail when its void volume fraction reaches 5%. Due to infinitesimal deformation of the solid half-sphere during the collision process, the material of the indenter is defined as rigid.
Vertical force versus displacement of the indenter curves from the experimental and numerical simulations is shown in Fig. 18. The structural deformation develops gradually with the increase of indentation depth. When the displacement of the indenter reaches approximately 0.16 m, the void volume fraction of the material near the symmetric centre of the upper plate is close to 5%, as shown in Fig. 15(a), and the vertical resistance of the side-shell structure reaches the maximum at this moment. With the continuous deformation, fracture occurs in the vicinity of the upper plate, which is shown by the deletion of elements with the void volume fraction reaching 5% in this numerical calculation, as shown in Fig. 15(b). The vertical force decreases dramatically. When the indenter penetration reaches approximately 0.31 m, the stiffeners attached to the lower plate and the lower plate itself withstand the load, and the vertical force increases gradually with the deformation development. The second peak of vertical force occurs at the moment when the lower plate breaks, and the distributions of void volume fraction before and after fracture are shown in Fig. 16. The total structural deformation at the termination time (excluding the reinforcing frame) is shown in Fig. 17. It can be seen that the error increases after the first peak point, which may be attributed to the increase of experimental uncertainty with the increasing number of structures involved in deformation after rupture occurs (Fig. 18). In addition, the following sources are not possible to include in the FE model, such as the residual stresses caused by welding, weld geometries, and roll direction of the material. Although the residual stresses were relaxed before the start of the experiment, it was done for the upper plate and not for the lower one (Table 3), which may cause the increasing error after the first peak point. In Hogström and Ringsberg (2012), the numerical results calculated using three different damage initiation (DI) criteria (shear, FLD, and FLSD) combined with a bilinear law for damage evolution (DE) are compared with experimental results, and the comparison of the best three combinations in various DI+DE models is shown in Fig. 19. Before the first peak, the calculated results of the GTN damage model are more consistent with the experimental results, as shown in Fig. 19. After that step, the numerical results calculated by all methods are slightly lower than the test results. In general, the numerical simulation results calculated by the GTN (a) Before fracture (b) After fracture Fig. 16 The distribution of void volume fraction on the lower plate before and after fracture

Conclusions
To better describe the intricate ductile failure process of NVA shipbuilding steel and accurately predict the structural deformation and damage during ship collision and grounding accidents, a new method to determine the critical void volume fraction is proposed based on mesoscopic damage model in this paper. A comparison between experiment and numerical simulations is used to validate the correctness of the method.
The following conclusions can be drawn: The GTN damage model is used in ship collision and grounding scenarios. Compared with the commonly used von Mises yield criterion, the mesoscopic damage model, as the foundation of the parameter determination method, links the material yield with damage, thereby reflecting the continuous deterioration of the material as the damage processes. Therefore, the GTN damage model is an ideal mathematical model describing the material rupture process. A new method for predicting the critical void volume fraction f c of the GTN damage model is proposed. The proposed critical void volume fraction identification method reduces the uncertainty of the GTN damage model parameter selection to some extent. A combination of the Voce and Swift equations is used to describe the stress-strain relationship of NVA mild steel to represent its ductile fracture process. The material parameters and the GTN damage parameters of NVA shipbuilding steel are determined by numerical comparison with a uniaxial tensile test, and reference experiment of a double-hull side-shell structure subjected to collision load is used to verify the correctness of the proposed method and the accuracy of the selected material parameters. The calculated results of the GTN damage model are also compared with the numerical results from Hogström and Ringsberg (2012), and the comparison results demonstrate the accuracy and applicability of the proposed method in this paper.
Combined with the connective strain hardening curve and the determinate material parameters, the GTN damage model with the new method will accurately describe the ductile fracture of NVA mild steel to better predict structural responses in collision accidents.
Open Access This article is distributed under the terms of the Creative Comm ons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Fig. 19 Comparison of FE analyses using GTN damage model and three failure criteria from Hogström and Ringsberg (2012)