A Modiﬁed Phase-Field Damage Model for Metal Plasticity at Finite Strains: Numerical Development and Experimental Validation

: Steel structures are designed to operate in an elastic domain, but sometimes plastic strains induce damage and fracture. Besides experimental investigation, a phase-field damage model (PFDM) emerged as a cutting-edge simulation technique for predicting damage evolution. In this paper, a von Mises metal plasticity model is modified and a coupling with PFDM is improved to simulate ductile behavior of metallic materials with or without constant stress plateau after yielding occurs. The proposed improvements are: (1) new coupling variable activated after the critical equivalent plastic strain is reached; (2) two-stage yield function consisting of perfect plasticity and extended Simo-type hardening functions. The uniaxial tension tests are conducted for verification purposes and identifying the material parameters. The staggered iterative scheme, multiplicative decomposition of the deformation gradient, and logarithmic natural strain measure are employed for the implementation into finite element method (FEM) software. The coupling is verified by the ‘one element’ example. The excellent qualitative and quantitative overlapping of the force-displacement response of experimental and simulation results is recorded. The practical significances of the proposed PFDM are a better insight into the simulation of damage evolution in steel structures, and an easy extension of existing the von Mises plasticity model coupled to damage phase-field.


Introduction
Engineering steel structures are designed to satisfy the demands of structural safety [1,2]. During their use, structures are intended to be exposed to predicted loading conditions depending on their purpose [3]. However, in some cases, due to unpredicted loading conditions (static, dynamic, or cyclic loading [4]), environments (corrosive [5,6] or high temperature [7]), or deviations in the design process, a non-permissible deformation and strain in the structure can be noticed [8]. Such behavior can lead to the initiation and evolution of damage, which often terminates in the structure's failure [9].
A cracking during the fracture phenomena in steel structures is delicate for numerical simulation [10][11][12]. Based on the Griffith theory, crack growth is related to a balance between bulk elastic energy and surface energy for brittle materials, which can be extended by plastic deformation energy and plastic dissipative energy for ductile materials [13,14]. A phasefield damage model (PFDM) can be used to overcome the limitations of the Griffith theory (a preexisting crack and a well-defined crack path) [15] and to show that a diffusive crack modeling is suitable for the numerical modeling of fracture [16][17][18]. Most PFDMs are based on variational principles [15]. Miehe et al. in [16,17] presented such a framework for diffusive fracture based on the phase-field approach and considered its numerical implementation by operator split scheme. A local history field is defined to govern the evolution of the crack field. Ambati et al. [19] discussed the problem of phase-field modeling of brittle fracture based on the Griffith theory. They proposed a hybrid formulation of a 'staggered' scheme. Molnar et al. [20] also proposed a phase-field model for brittle fracture based on the rate-independent variational principle of diffuse fracture. The displacement field and phase-field are solved separately by the 'staggered' approach. The phase-field variable separates the damaged, and virgin material states smoothly [21]. It is a scalar value, zero for virgin, and one for damaged material [20]. The displacement and phase field can be coupled in the 'monolithic' or 'staggered' approach [22]. The 'staggered' algorithm needs a well-defined stopping criterion [19].
Various authors developed and implemented PFDM in commercial or in-house finite element method (FEM) software. Azinpour et al. [23] considered the analogy between temperature and failure models and proposed unified and straightforward implementation of PFDM into Abaqus FEM software. Liu et al. [22] used commercial FEM software Abaqus to explore the monolithic and staggered coupling approaches for the phase-field model. They used user material and user element subroutines. The coupling between plasticity and PFDM is essential for the efficient modeling of cracking processes [13]. Plasticity is related to the development of inelastic strains, while damage is associated with reducing the material's stiffness [24]. A contribution to the energy dissipation is related to plastic strains [13]. Ambati et al. in [18] presented a phase-field model for ductile fracture capable of capturing the behavior of J2-plasticity material and crack initiation, propagation, and failure. Ambati et al. in [25] extended the phase-field model to the finite strain regime. Fracture is controlled by a critical scalar value of plastic strain. This model allowed simulation of various phenomena such as necking and fracture of flat specimens. They used the S355 type of steel for validation purposes. Badnava et al. [26] proposed a phase-field and ratedependent plasticity coupling by energy function. The influence of plastic strain energy on crack propagation was defined by a threshold variable. Miehe et al. [27] presented a variational formulation for the phase-field modeling of ductile fracture for large strains. They introduced independent length scales for the plastic zone and cracks to guarantee mesh objectivity in post-critical ranges. Paneda et al. in [28] presented the possibility of the phasefield approach for fracture as a suitable solution for hydrogen-assisted cracking to predict catastrophic failures in corrosive environments. A fine mesh is necessary for the smooth phase-field distribution what can be computationally expensive [21]. Zhang et al. [29] investigated the accuracy of crack path simulation by the small length scale, leading to an unrealistic force-displacement relationship. Seles et al. [12,21] investigated and presented a stopping criterion based on the residual norm's control within a fracture analysis staggered scheme. Ribeiroa et al. [30] presented possibilities of the Abaqus damage plasticity model for simulation of S355NR+J2 experimentally investigated specimens. They obtained the specific behavior of this type of steel in the yielding zone [31].
The main research findings presented in this article are enhanced simulation of steel structure behavior by the coupled PFDM and von Mises plasticity model for ductile fracture presented by Ambati et al. [18,25], and Miehe et al. [16,17,32]. The improvements of coupled multifield three-dimensional finite element and Simo's hardening function for plasticity are implemented into the in-house FEM software PAK developed at the Faculty of Engineering, University of Kragujevac, Serbia. In Section 2.1, the overview of governing equations evolution for the PFDM and the von Mises plasticity with Simo hardening function is presented. The main improvements are: • the modified coupling variable presented in Section 2.2, which considers the damage influence induced by plastic strains after the saturation hardening stress is achieved, and • the two-stage yield function presented in Section 2.3 for the simulation of metallic materials behavior, which exhibits the stress plateau after yielding occurs. It consists of (1) perfect plasticity and (2) extended Simo-type yield hardening function. The evolution of the extended Simo-type yield function from its basic form for simple implementation into the standard von Mises plasticity constitutive model is given.
In Section 2.4, the large-strain theory based on the multiplicative decomposition of the deformation gradient, and logarithmic natural strain measure is given in the algorithmic form. The theory is adopted to fit the staggered iterative scheme for the displacement and damage phase-field solution. In Section 2.5, the FEM implementation details including, the staggered coupling scheme, are presented. Section 3 is focused on the verification of staggered iterative scheme by benchmark one element example available in the literature and the experimental investigation of S355J2+N steel flat dog-bone specimens and validation of PFDM and von Mises plasticity model. The material parameters are identified by the calibration of the force-displacement diagram comparing the experimental and numerical results. By comparing the force-displacement response of the experimental investigation and simulation results presented an excellent qualitative and quantitative prediction, while the equivalent plastic strain field compared to the deformed configuration of the specimen gives good qualitative comparison results. At the end of the article, in Section 4, the main conclusions are presented.

Overview of Governing Equations Evolution
In this subsection, the governing equations for coupled damage phase-field and plasticity are derived according to the literature [15][16][17]20,28,[32][33][34][35] to clarify the proposed modifications. According to Griffith's theory, a fracture is defined by a criterion based on the equilibrium of the surface energy and the elastic energy stored in the material. It is possible to predict a crack initiation for existing cracks, but the nucleation and the crack propagation is not possible. There are two primary methods for modeling crack propagation in structures: (a) discrete and (b) diffuse [20]. The diffuse method considers crack as smeared damage. Francfort and Marigo [15] proposed a variational fracture model based on minimizing an energy functional for the displacement field and discontinuous crack set. The model of Bourdin et al. [33] is regularized with a smeared crack by the introduction of a phase-field to describe fully broken and intact material phases.
To introduce the phase-field model as a type of diffuse method for crack modeling, let us consider a bar given in Figure 1 with a constant cross-section. A damage phase-field variable d along the coordinate x of the bar can be formulated as local discontinuity for sharp crack topology as follows [16,17,20,28] however, for the diffusive crack topology, the damage can be given as an exponential function of the bar length in the form where l c is defined as a characteristic length-scale parameter. This formulation converges to Equation (1), when l c → 0 . By following Miehe et al. [17] formulation for cracks in one-dimensional solids, and the extension of the regularized crack functional to multi-dimensional problems as  [17,20]. By following Miehe et al. [17] formulation for cracks in one-dimensional solids, and the extension of the regularized crack functional to multi-dimensional problems as a crack surface density function γ per unit volume is defined as [16,17,20,28,32] where ∇ is the gradient operator. For the application of the phase-field model to the ductile behavior of the materials, the internal potential energy density ψ for the ductile fracture is considered as the sum of elastic ψ E (ε E , d) and plastic ψ P (ε P , d) energy density, fracture surface energy density ϕ S (d) and plastic dissipated energy density ϕ P (ε P , d) as in [27] where ε E is the elastic strain tensor, and ε P is the equivalent plastic strain. Let us define each term in Equation (5). The elastic energy density of virgin material ψ E 0 is multiplied by degradation function g(d) to define the elastic energy density ψ E as [20,27,28] where C 0 is the fourth-order elastic constitutive matrix, and σ 0 is the Cauchy stress tensor of an undamaged solid. Similarly, the "damaged" Cauchy stress σ is given in the following form [17] By using Equation (3), the fracture surface energy Φ S at the crack surface S is defined as [16,17] where the fracture surface energy density dissipated by the formation of the crack is defined as In Equations (8) and (9), the Griffith-type critical fracture energy release rate G c is used, also known as the fracture toughness of the material described as the amount of energy required to produce a unit area of fracture surface [33,34]. The plastic energy density for Simo hardening is [17] where σ yv is the initial yield stress, σ y0,∞ is the saturation hardening stress, n is the hardening exponent, and H is the hardening modulus [36]. The plastic dissipated energy density is [35] ϕ P (ε P , d) = g(d)σ yv ε P .
By using the Equations (5), (6) and (9)-(11) the total internal potential energy Ψ functional is defined as [35] where the authors introduced a critical fracture energy release rate per unit volume as G V = G c /l c . The variation of the internal potential energy in Equation (12) over the elastic strain, damage and equivalent plastic strain is given as [35] and for the Simo hardening function [17], one can obtain [27,35] where ε P is the plastic strain tensor, and g (d) is the derivative of the degradation function g(d) over d. A variation of the external potential energy W ext is known as [37] where b is a body force field per unit volume, h is a boundary traction per unit area, and u is the displacements vector. The equilibrium of the internal (14) and external (15) potential energy for the Simo hardening function gives [35] V σ : By the application of total derivatives of the following terms and by using the Gauss theorem, it can be obtained [27,35] where n is the unit outer normal to the surface A, and the internal potential energy density given in Equation (5) is The Neumann-type boundary conditions are [20] σ · n − h = 0, (23) ∇d · n = 0, what leads to the governing balance equations of the coupled PFDM-von Mises plasticity problem for Simo's hardening function [35] ∇d · n = 0, where the equivalent stress is defined as σ eq = σ 0 : ∂ε P ∂ε P .

Modifications of Coupling Variable
The degradation function and its derivative over d are proposed by Ambati et al. [18] for the phase-field damage modeling of ductile fracture as where p is the coupling variable. The same degradation function in Equation (28) can be used for the brittle fracture by setting p = 1 what will be used for the one element example in Section 3. The coupling variable can be defined to depend on the critical value of equivalent plastic strain ε crit P [18]. In this paper, the authors propose the modification of the Ambati et al. [18] because the material is considered to be intact (undamaged) until the equivalent plastic strain achieves the critical value ε P = ε crit P . The critical value of the equivalent plastic strain can be registered in the experimental stress-strain diagram, when the stress starts to decrease ( Figure 2).

Modifications of Coupling Variable
The degradation function and its derivative over d are proposed by Ambati et a [18] for the phase-field damage modeling of ductile fracture as ( ) ( ) where p is the coupling variable. The same degradation function in Equation (28) can b used for the brittle fracture by setting 1 p = what will be used for the one element exam ple in Section 3. The coupling variable can be defined to depend on the critical value o equivalent plastic strain crit P ε [18]. In this paper, the authors propose the modification o the Ambati et al. [18] because the material is considered to be intact (undamaged) unt the equivalent plastic strain achieves the critical value  At that critical point C (Figure 3b), the material can be considered damaged due t the plastic strains and the damage-plasticity coupling variable is activated. Therefore, th coupling variable p can be defined for equivalent plastic strain as it is given in Figure   Figure 2. New coupling variable p (continuous line) in relation to the equivalent plastic strain ε P compared to Ambati et al. [18].
At that critical point C (Figure 3b), the material can be considered damaged due to the plastic strains and the damage-plasticity coupling variable is activated. Therefore, the coupling variable p can be defined for equivalent plastic strain as it is given in Figure 3 by the function Metals 2021, 11, x FOR PEER REVIEW 8 of 27 The stored internal potential energy density ψ is considered as the elastic energy density 0 E ψ because the influence of the plastic part is taken into account by coupling variable p , so [25] where t ψ is the previously stored internal potential energy density.

Two-Stage Yield Hardening Function
For the simulation of the metallic material's behavior, which exhibits stress plateau after yielding occurs, such as S355J2+N steel, the extended two-stages yielding function is necessary to describe the idealized response given in Figure 3 (continuous line). In the first stage, the yielding occurs, the plastic strain increases, while the stress is constant. In the second stage, the stress increases nonlinearly until the saturation hardening stress 0, y σ ∞ . After the stress achieves the maximal value at the end of the second stage, the stress decreases due to the phase-field damage model influence on the material (dashed line). To simulate the described behavior, 'perfect plasticity' is employed for the first stage of loading ( 0 P P ε ε < ) until the plastic strain 0 P ε is achieved. In the first stage, the yield condition is given in the following form [37] 0 while in the second period, the yield condition is defined based on Equation (28) as The stored internal potential energy density ψ is considered as the elastic energy density ψ E 0 because the influence of the plastic part is taken into account by coupling variable p, so [25] where t ψ is the previously stored internal potential energy density.

Two-Stage Yield Hardening Function
For the simulation of the metallic material's behavior, which exhibits stress plateau after yielding occurs, such as S355J2+N steel, the extended two-stages yielding function is necessary to describe the idealized response given in Figure 3 (continuous line). In the first stage, the yielding occurs, the plastic strain increases, while the stress is constant. In the second stage, the stress increases nonlinearly until the saturation hardening stress σ y0,∞ . After the stress achieves the maximal value at the end of the second stage, the stress decreases due to the phase-field damage model influence on the material (dashed line). To simulate the described behavior, 'perfect plasticity' is employed for the first stage of loading (ε P < ε P0 ) until the plastic strain ε P0 is achieved. In the first stage, the yield condition is given in the following form [37]  while in the second period, the yield condition is defined based on Equation (28) as where the yield stress equal to equivalent stress in the extended Simo-type yield function is given as The complete two-stage yield function shown in Figure 3 by continuous line can be defined by the equation If the yield function (35) is less than zero, the solution is elastic. If the condition is violated, the equivalent plastic strain increment ∆ε P of the function f y = 0 must be determined in an iterative Newton-Raphson procedure given in the stress integration algorithm in the next subsection. After the correct value of equivalent plastic strain ε P is computed, the deviatoric stress tensor, the total stress tensor, the elastic strain tensor, and the elastic-plastic constitutive matrix must be updated.
The yield stress term given by Equation (34) in the extended Simo-type yield function is derived from its basic form obtained from in Equation (27) as The Equation (34) can be transformed into the form To provide a simple implementation into the existing von Mises constitutive model for metal plasticity, the basic form of Simo's yield function σ yS in Equation (36) can be transformed by subtraction of the stress increment ∆σ yS for the plastic strain ε P0 , when the constant stress was registered during yielding given as To obtain the stress function in Equation (37) it is necessary to calculate the difference between the basic form of Simo's hardening function in Equation (36) and the stress value given by Equation (38) as which finally gives By comparing the yield stress in Equations (37) and (40), it can be noticed that the second term in Equation (40) is multiplied by e −nε P0 , so to propose the same form as Equation (37), and to allow the possibility of using the Equation (39) for the implementation purpose, the coefficient in the second term in Equation (40), σ y0,∞ − σ yv needs to be multiplied by e nε P0 what will give the equivalent equation to Equation (37).

Stress Integration Algorithm for von Mises Large Strain Plasticity
In this section, an overview of the well-known stress integration algorithm for von Mises large strain plasticity is presented with the two new steps 11 and 12 which are necessary for solving the PFDM governing equation given in Equation (26). The value of elastic strain energy ψ E 0 is calculated at the integration point level as well as the coupling variable p. In the following equation, the complete algorithm is given to provide the implementation. The deviatoric strain can be obtained using the multiplicative decomposition of the deformation gradient [37] where F E and F P are the elastic and plastic deformation gradient, respectively. The elastic deformation gradient can also be decomposed into the isochoric deformation gradient F E and volumetric portion as [37] so, the elastic left Cauchy-Green strain tensor b E can be calculated as [37,38] The elastic deviatoric strain e E can be calculated using the Hencky strain measure h E as [37,38] and the mean strain e m , in that case, is [37,38] The total stress tensor can be decomposed on the elastic deviatoric S E and the volumetric part σ m as [37] σ = S E + σ m I, where I is the unit tensor. The elastic deviatoric stress can be defined as [37] and the mean stress is [37] where shear and bulk modulus are while E is the Young's modulus, and ν is the Poisson's ratio. The detailed algorithm of the von Mises plasticity for large strain problems, is given in details below [37]: t-time at the beginning of time step; ∆t-time increment Initial conditions (save at the integration point level): 3.

5.
Check for yielding: If the condition is satisfied, the solution is S E = S * E and ∆ε P = 0, and go to 7. 6.
Find equivalent plastic strain increment ∆ε P of the function f y (∆ε P ) = 0 7.
Update of left Cauchy-Green strain tensor: 8.
Mean stress and total stress: 9.

Implementation into FEM Software
In this section, the discretization using standard Lagrange finite elements is considered and the standard Galerkin finite element method is used. The nodal displacements and phase-field damage variable are unknowns that need to be determined. These finite elements are also known as "multifield" finite elements extensively applied in multiphysics FE simulations.

Finite Element Discretization
The interpolation matrices for displacement N u and damage N d and derivatives B d and B u are given as follows [20]: The damage phase field value at an integration point is described as [20] where d is the damage phase-field vector of nodal values. The local damage gradient which can also be referred to as "damage strain" is [20] The total strain vector ε is interpolated in terms of the nodal displacements u [20,37] The internal f int e and external forces f ext e are [20,37] f int The residue vector r d e for the phase-field degrees is [20,37] The tangent matrices for damage K d e and displacement K u e field are [20,37] The linear equation system can be solved by using of Newton-Raphson algorithm [20]

Staggered Solution Strategy
By following Miehe et al. [16], the weak formations of mechanical field given in Equations (25) and (27) and phase-field given in Equation (26) have to be solved by the staggered algorithm shown in Figure 4. The Equation (27) is the Simo yield hardening function used as a part of the local stress integration algorithm for the computation of the mechanical field. At the beginning of the iterative procedure (t 0 ), both, displacement and damage vectors are equal to the vectors from the previous time step. Initial conditions at the structure level are u (0) = t u; d (0) = t d.
Metals 2021, 11, x FOR PEER REVIEW

Staggered Solution Strategy
By following Miehe et al. [16], the weak formations of mechanical field given tions (25) and (27) and phase-field given in Equation (26) have to be solved by t gered algorithm shown in Figure 4. The Equation (27)   Both fields are solved simultaneously by minimization of the two residual eq for the displacement u e r and the damage phase-field d e r to determine the vecto next time step [18,20] The staggered iterative scheme across the time steps is given in Figure 4, whe be observed that the procedure starts with both displacement and damage vecto to zero. The Newton-Raphson iterative procedure is given with the convergence and the implementation details below [20,37]: Both fields are solved simultaneously by minimization of the two residual equations for the displacement r u e and the damage phase-field r d e to determine the vectors in the next time step [18,20] The staggered iterative scheme across the time steps is given in Figure 4, where it can be observed that the procedure starts with both displacement and damage vectors equal to zero. The Newton-Raphson iterative procedure is given with the convergence criterion and the implementation details below [20,37]: Stress integration σ Stiffness matrices Displacement and damage increment, update of displacement and damage If convergence criteria are not satisfied, go to step B.
If r u(i) ≤ tol and r d(i) ≤ tol go to next time step A.

One Element-Brittle Fracture Benchmark Example
The simplest model that can be used to verify and understand the proposed staggered iterative scheme (Section 2.5.2) for coupling of damage phase field and displacement field is suggested by Molnar and Gravouil in [20] for the brittle fracture. The model is one three-dimensional hexahedral element of unit dimensions (1 × 1 × 1 mm).
The parameters necessary for brittle fracture simulation are the same as in [9]: the Young's modulus E = 210 GPa, the Poisson's ratio υ = 0.3, the critical energy release rate G V = 5 · 10 −2 GPa and the length scale parameter l c = 0.1 mm. The plasticity material parameters are large enough so the material is in the elastic regime, while the coupling variable is p = 1 to simulate the brittle fracture [9]. The analytical stress σ vs. strain ε, as well as damage d vs. strain ε relationships, are given as follows [9]: where c 22 is the element of the elastic constitutive matrix C 0 . The nodes on the bottom of the cube are constrained in all three directions, while the top nodes are allowed to slide vertically. The loading is realized in a displacement control regime in 1000 time steps until the total displacement of 0.1 mm. The obtained results are quantitatively compared to the analytical results [9]. The comparison given in Figures 5 and 6 confirms the functionality of the proposed iterative scheme and its correctness.

Experimental Investigation and FEM Simulation of S355J2+N Specimens
The S355J2+N is widely used in engineering structures due to the good weldability and machinability. It also exhibits constant stress plateau after yielding occurs, so it is chosen as a representative steel type to validate the two-stage hardening yield function. This material is also used in experimental investigation and simulation of fracture by various authors [25,30,39]. For the validation purposes of modified PFDM, three steel S355J2+N specimens are investigated by using servo-hydraulic testing machine-EHF-

Experimental Investigation and FEM Simulation of S355J2+N Specimens
The S355J2+N is widely used in engineering structures due to the good weldability and machinability. It also exhibits constant stress plateau after yielding occurs, so it is chosen as a representative steel type to validate the two-stage hardening yield function. This material is also used in experimental investigation and simulation of fracture by various authors [25,30,39]. For the validation purposes of modified PFDM, three steel

Experimental Investigation and FEM Simulation of S355J2+N Specimens
The S355J2+N is widely used in engineering structures due to the good weldability and machinability. It also exhibits constant stress plateau after yielding occurs, so it is chosen as a representative steel type to validate the two-stage hardening yield function. This material is also used in experimental investigation and simulation of fracture by various authors [25,30,39]. For the validation purposes of modified PFDM, three steel S355J2+N specimens are investigated by using servo-hydraulic testing machine-EHF-EV101 K3-070-0A (Shimadzu Corporation, Tokyo, Japan), with force ±100 kN and stroke ±100 mm. The specimen's chemical composition given in Table 1.  [40] and ASTM E8M-01 [41] at room temperature (23 ± 5 • C) in constant stroke control rate of 4 mm/min. Specimen's shape and dimensions are given in Figure 7. For the elongation measurement and identification of Young modulus, the extensometer MFA25 (MF Mess-& Feinwerktechnik GmbH, Velbert, Germany), with a gauge length of 50 mm is used as given in Figure 8.   [40] and ASTM E8М-01 [41] at room temperature ( 23 5 C ±°) in constant stroke control rate of 4 mm/min. Specimen's shape and dimensions are given in Figure 7. For the elongation measurement and identification of Young modulus, the extensometer MFA25 (MF Mess-& Feinwerktechnik GmbH, Velbert, Germany), with a gauge length of 50 mm is used as given in Figure 8.  The investigated specimens are presented in Figure 9 after the experiment. The forcedisplacement responses are recorded and the comparison of the obtained results is given in Figure 10. The response of "Specimen 3" is selected as the representative for the PFDM validation purpose.   [40] and ASTM E8М-01 [41] at room temperature ( 23 5 C ±°) in constant stroke c trol rate of 4 mm/min. Specimen's shape and dimensions are given in Figure 7. For elongation measurement and identification of Young modulus, the extensometer MFA (MF Mess-& Feinwerktechnik GmbH, Velbert, Germany), with a gauge length of 50 m is used as given in Figure 8.  The investigated specimens are presented in Figure 9 after the experiment. The for displacement responses are recorded and the comparison of the obtained results is giv in Figure 10. The response of "Specimen 3" is selected as the representative for the PFD validation purpose. The investigated specimens are presented in Figure 9 after the experiment. The forcedisplacement responses are recorded and the comparison of the obtained results is given in Figure 10. The response of "Specimen 3" is selected as the representative for the PFDM validation purpose.  The FE model is prepared for the straight part of the specimen same as the gauge length (50 mm), according to the specimen's dimensions. One-eighth of the specimen is modeled due to the existence of three symmetry planes. Dimensions of the FE model are 25 × 6.25 × 2.5 mm. The geometrical imperfection necessary to trigger the plastic deformation process in a zone of 10 mm (L2) from the middle of the specimen, where necking is expected, is prescribed as 0.01% a linear decrease of the specimen width D and thickness ( Figure 11). The FE model is created using 2100 standard full integrated 8-node hexahedral finite elements with mesh refinement in the expected necking zone. Boundary conditions include the constraint of nodes in symmetry planes in a direction perpendicular to the symmetry plane they belong to. The FE model tensile loading is applied to the top surface nodes by displacement increment of 0.02 mm for 350 steps.  The FE model is prepared for the straight part of the specimen same as the gau length (50 mm), according to the specimen's dimensions. One-eighth of the specimen modeled due to the existence of three symmetry planes. Dimensions of the FE model 25 × 6.25 × 2.5 mm. The geometrical imperfection necessary to trigger the plastic def mation process in a zone of 10 mm (L2) from the middle of the specimen, where necki is expected, is prescribed as 0.01% a linear decrease of the specimen width D and thi ness ( Figure 11). The FE model is created using 2100 standard full integrated 8-node h ahedral finite elements with mesh refinement in the expected necking zone. Bounda conditions include the constraint of nodes in symmetry planes in a direction perpendi lar to the symmetry plane they belong to. The FE model tensile loading is applied to top surface nodes by displacement increment of 0.02 mm for 350 steps. The FE model is prepared for the straight part of the specimen same as the gauge length (50 mm), according to the specimen's dimensions. One-eighth of the specimen is modeled due to the existence of three symmetry planes. Dimensions of the FE model are 25 × 6.25 × 2.5 mm. The geometrical imperfection necessary to trigger the plastic deformation process in a zone of 10 mm (L 2 ) from the middle of the specimen, where necking is expected, is prescribed as 0.01% a linear decrease of the specimen width D and thickness ( Figure 11). The FE model is created using 2100 standard full integrated 8-node hexahedral finite elements with mesh refinement in the expected necking zone. Boundary conditions include the constraint of nodes in symmetry planes in a direction perpendicular to the symmetry plane they belong to. The FE model tensile loading is applied to the top surface nodes by displacement increment of 0.02 mm for 350 steps. Metals 2021, 11, x FOR PEER REVIEW 19 of 27 Figure 11. Finite element mesh, imperfection, loading, and boundary conditions.
The material parameters used for simulations are given in Table 2:   Figure 12 shows the deformed S355J2+N specimen after the experimental investigation along with the equivalent plastic strain field obtained from the PFDM FEM simulation. The distribution of the equivalent plastic strain field qualitatively simulates the deformed configuration of the experimentally obtained deformations. The equivalent plastic strain field obtained by phase-field plasticity coupled simulation is localized in a zone where the fracture occurs. The strong relationship between the damage field and the equivalent plastic strain field obtained by the coupled phase-field simulation given in Figure 13 imposes that damage is a governing phenomenon that leads to the fracture of the specimen. The equivalent plastic strain field for the pure plasticity without phase-field (Figure 13a) and the PFDM simulation (Figure 13b) are presented to show the influence of the damage field on the The material parameters used for simulations are given in Table 2 Figure 12 shows the deformed S355J2+N specimen after the experimental investigation along with the equivalent plastic strain field obtained from the PFDM FEM simulation. The distribution of the equivalent plastic strain field qualitatively simulates the deformed configuration of the experimentally obtained deformations. The equivalent plastic strain field obtained by phase-field plasticity coupled simulation is localized in a zone where the fracture occurs. The material parameters used for simulations are given in Table 2:   Figure 12 shows the deformed S355J2+N specimen after the experimental investigation along with the equivalent plastic strain field obtained from the PFDM FEM simulation. The distribution of the equivalent plastic strain field qualitatively simulates the deformed configuration of the experimentally obtained deformations. The equivalent plastic strain field obtained by phase-field plasticity coupled simulation is localized in a zone where the fracture occurs. The strong relationship between the damage field and the equivalent plastic strain field obtained by the coupled phase-field simulation given in Figure 13 imposes that damage is a governing phenomenon that leads to the fracture of the specimen. The equivalent plastic strain field for the pure plasticity without phase-field (Figure 13a) and the PFDM simulation (Figure 13b) are presented to show the influence of the damage field on the The strong relationship between the damage field and the equivalent plastic strain field obtained by the coupled phase-field simulation given in Figure 13 imposes that damage is a governing phenomenon that leads to the fracture of the specimen. The equivalent plastic strain field for the pure plasticity without phase-field (Figure 13a) and the PFDM simulation (Figure 13b) are presented to show the influence of the damage field on the localization of the plastic strains. In Figure 13a, the plastic strains are localized in the middle of the model, with the minimal difference between the minimal and maximal value, which does not suggest the fracture zone location. On the other side, the distribution of the damage field given in Figure 13c corresponds to the equivalent plastic strain field in Figure 13b, so it can be considered as a generator of the fracture process. localization of the plastic strains. In Figure 13a, the plastic strains are localized in the middle of the model, with the minimal difference between the minimal and maximal value, which does not suggest the fracture zone location. On the other side, the distribution of the damage field given in Figure 13c corresponds to the equivalent plastic strain field in Figure 13b, so it can be considered as a generator of the fracture process. The comparison of the force-displacement relationship between experimental and simulation results is given in Figure 14. By comparing the diagrams, one can notice that 'pure' plasticity without damage cannot follow the behavior recorded by the experiment after the maximal force is achieved. Using the PFDM approach and the proposed modified coupling variable p , the influence of plastic strain development is activated after the loading force attains the maximum value and starts to decrease. The coupling variable p linearly increases as described in Figure 3 and simultaneously, the force starts to decrease until the specimen's fracture. In Figure 14, the value of the damage is given in the middle of the specimen in relation to the displacement of the specimen. It can be noticed that the damage is zero until the critical value of plastic strain is achieved. After the damage starts to increase, the force decreases following the slope of the damage change. The comparison of the force-displacement relationship between experimental and simulation results is given in Figure 14. By comparing the diagrams, one can notice that 'pure' plasticity without damage cannot follow the behavior recorded by the experiment after the maximal force is achieved. Using the PFDM approach and the proposed modified coupling variable p, the influence of plastic strain development is activated after the loading force attains the maximum value and starts to decrease. The coupling variable p linearly increases as described in Figure 3 and simultaneously, the force starts to decrease until the specimen's fracture. In Figure 14, the value of the damage is given in the middle of the specimen in relation to the displacement of the specimen. It can be noticed that the damage is zero until the critical value of plastic strain is achieved. After the damage starts to increase, the force decreases following the slope of the damage change.
The element's dimension of the coarse FE mesh along the specimen length is 0.5 mm. This mesh has been used in previous simulations to show that the force-displacement response can be obtained even for the coarse FE meshes. However, if we want to simulate the evolution of damage phase-field with moving interface, the FE mesh must be refined further in the zone where the crack is expected. The dimension of the first two rows of elements (1.0 mm of specimen length) is reduced 4 times for the medium mesh (the element length-0.125 mm), while for the fine mesh the reduction is 10 times (the element length-0.05 mm). The force-displacement responses for the same material parameters given in Table 2 are presented in Figure 15, where it can be noticed that the finer mesh gives softer response in post-critical zone as suggested in Ambati et al. [18]. The evolution of damage field for the post-critical behavior is given in Figure 16 as well as the equivalent plastic strain field development in Figure 17. As it can be noticed, both the damage field and the equivalent plastic strain field evolves in the cracking zone of the specimen. For the visualization purpose, the further research will be criteria of "element death" which need to be satisfied to remove the elements as suggested in [19]. The element's dimension of the coarse FE mesh along the specimen length is 0.5 mm. This mesh has been used in previous simulations to show that the force-displacement response can be obtained even for the coarse FE meshes. However, if we want to simulate the evolution of damage phase-field with moving interface, the FE mesh must be refined further in the zone where the crack is expected. The dimension of the first two rows of elements (1.0 mm of specimen length) is reduced 4 times for the medium mesh (the element length-0.125 mm), while for the fine mesh the reduction is 10 times (the element length-0.05 mm). The force-displacement responses for the same material parameters given in Table 2 are presented in Figure 15, where it can be noticed that the finer mesh gives softer response in post-critical zone as suggested in Ambati et al. [18]. The evolution of damage field for the post-critical behavior is given in Figure 16 as well as the equivalent plastic strain field development in Figure 17. As it can be noticed, both the damage field and the equivalent plastic strain field evolves in the cracking zone of the specimen. For the visualization purpose, the further research will be criteria of "element death" which need to be satisfied to remove the elements as suggested in [19].          Finally, the von Mises model with the same hardening function, but without coupled damage field (PFDM), cannot capture the post-critical behavior resulting from the material deterioration. The S355J2+N type of steel exhibits the specific plateau after the yielding occurs and it is captured by the proposed two-stage hardening function. The standard Simo's hardening function cannot be used for the simulation of such materials, as shown in Figure 18. The standard Simo's hardening function cannot follow the experimental forcedisplacement response at the same quantitatively and qualitatively level as the proposed two-stage hardening yield function. damage field (PFDM), cannot capture the post-critical behavior resulting from the m rial deterioration. The S355J2+N type of steel exhibits the specific plateau after the yie occurs and it is captured by the proposed two-stage hardening function. The stan Simo's hardening function cannot be used for the simulation of such materials, as sh in Figure 18. The standard Simo's hardening function cannot follow the experim force-displacement response at the same quantitatively and qualitatively level as the posed two-stage hardening yield function. Figure 18. Comparison of the standard and the proposed two-stage hardening function for the same material parameters given in Table 1.

Conclusions
(1) The purpose of this article is to offer a better insight into the damage of material fracture of structures, propose necessary modifications of the existing PFDM allow better control over the damage simulation process. Simo-type hardening) to offer more realistic simulations of metallic materials b ior, which exhibit constant stress plateau after yielding. (4) To control the start of the damage field development, the modification of the pling between the plastic strain and the damage field is determined by the cou variable, which starts to increase linearly after the equivalent plastic strain ach the critical value. (5) The main differences and advantages of the proposed method are possibilities to trol: (a) the onset of hardening after the initial constant stress plateau is ende the initiation of damage phase-field development due to the plastic strain dev ment, (c) the distribution of critical fracture energy release rate in damaged zon Figure 18. Comparison of the standard and the proposed two-stage hardening function for the same material parameters given in Table 1.

Conclusions
(1) The purpose of this article is to offer a better insight into the damage of materials and fracture of structures, propose necessary modifications of the existing PFDM, and allow better control over the damage simulation process. Simo-type hardening) to offer more realistic simulations of metallic materials behavior, which exhibit constant stress plateau after yielding. (4) To control the start of the damage field development, the modification of the coupling between the plastic strain and the damage field is determined by the coupling variable, which starts to increase linearly after the equivalent plastic strain achieves the critical value. (5) The main differences and advantages of the proposed method are possibilities to control: (a) the onset of hardening after the initial constant stress plateau is ended, (b) the initiation of damage phase-field development due to the plastic strain development, (c) the distribution of critical fracture energy release rate in damaged zone. (6) The successful implementation of the staggered coupling scheme has been verified by one element benchmark example from literature. The same results have been obtained for both stress-strain response, but also the damage-strain relationship. (7) The main application of the implementation is shown by simulation of S355J2+N test specimen behavior investigated by the universal testing machine. The experimentally captured force-displacement response is compared to the FEM simulation by proposed modified PFDM and excellent results are achieved. Also, the evolutions of the equivalent plastic strain field and the damage phase field are presented and the zone of maximal plastic strain qualitatively corresponds to the main deformed zone of the experimentally investigated specimens.