Numerical Formulation of Anisotropic Elastoplastic Behavior Coupled with Damage Model in Forming Processes

Abstract

The present paper proposes a mathematical development of the plasticity and damage approaches to simulate sheet metal forming processes. It focuses on the numerical prediction of the deformation of the sheet metal during the deep drawing process when a crack appears. Anisotropic plasticity constitutive equations are proposed. A fully implicit integration of the coupling constitutive equations is used and leads to two nonlinear local scalar equations that are solved by Newton’s method. The developed model allows predicting the onset of cracks in sheet metals during cold forming operations. The numerical model is implemented in ABAQUS software using user-defined subroutines, which are VUMAT and UMAT. The accuracy of the anisotropic elastoplastic model fully coupled with ductile damage is evaluated using numerical examples.

1. Introduction

Cold forming processes are widely used in industrial engineering. They allow for the production of principally sheet workpieces of various domains such as automotive and aeronautic parts. These processes are generally characterized by a high degree of permanent deformation of the blank under the action of the punch and die. Many approaches have been developed to predict the mechanical behavior of sheet materials [1,2]. In addition, several numerical approaches have been developed and applied in the simulations of sheet metal forming processes such as deep drawing and stamping operations. The study of any sheet forming process requires a knowledge of mechanical behavior evolution until the damage of the blank [3,4]—the large irreversible deformation leading to high strain localization zones and then cracks due to the ductile damage. Indeed, the sheets present anisotropy, which must be considered in the simulations of sheet shaping [5,6]. However, one of the most used plasticity models is the classical Hill 1948 yield criterion. This criterion represents a generalized plasticity yield function for anisotropic materials when compared to J2 plasticity. In sheet metal forming, the workpiece tends to have discontinuities due to the large plastic deformations causing nucleation and the void’s growth. Then, cracks propagate until failure.

A wide variety of models for the failure mechanism is discussed in the literature [7,8]. They are divided principally into two types. The first type of models are developed in the context of fracture mechanics [9]. The second type are defined by Continuum Damage Mechanics (CDM) formulation, in which the damage field is characterized by an internal variable [10]. Recently, models incorporating damage and anisotropic yield functions have attracted the attention of many authors due to the high demand of accurate numerical models in metal forming industries. Damage models are mainly based either on the Gurson approach [11], or on CDM [10,12], among others. These approaches were developed to accurately simulate forming processes and predict damage in sheet workpieces. For the uncoupled approach, the damage field is calculated without considering its effect on the mechanical behavior model. However, in the coupled approach, the damage affects the overall constitutive equations of mechanical fields.

The Gurson model [11] is based on the localization of deformation in a specific necking zone. The constitutive relation related to the Gurson theory consists of micro structural criterion in porous materials. The variation of the void’s volume during failure evolution is considered and added to the yield function expression by using the void density parameter. In the same context, the Gurson, Needleman, and Tvergaard (GNT) model [13,14] takes into account the effect of the void’s coalescence on the acceleration of the failure process. In the literature, many researchers are interested in the analysis of the growth and coalescence of voids during ductile damage [15,16,17,18].

Furthermore, CDM describes the evolution of the damage variable from a continuum mechanics approach. The implementation of the CDM model combined with the anisotropic elastoplastic constitutive equations require the use of advanced resolution schemes. Moreover, the elastoplastic models fully coupled with CDM include the effects of isotropic/kinematic hardening. The nonlinear isotropic/kinematic hardening and the consistent elastoplastic tangent modulus were acquired by [19]. In this work, the nonlinear material behavior is the result of the plasticity and the ductile damage phenomena. Indeed, numerical integration leads to four equations, which are two scalars and two tensors. De Souza et al. [20,21] used the finite strain extensions of this model. Nonlinear kinematic hardening is taken into account; however, the resulting return mapping schemes demand the same number of equations [22]. The backward Euler scheme is used to integrate the rate of constitutive relations. The consistent tangent operator was derived using the exact linearization of the algorithm.

In the literature, many formulations of the damage dissipation potential are proposed For ductile damage [10,12]. The continuum ductile damage model of Lemaitre [23] was usually used [19,22,24]. In the same context, a damage model was proposed [25] and used as a function of the accumulated plastic strain. In this model, a general nonlinear law of damage evolution is applied. In addition, Badreddine et al. [26] and Saanouni et al. [27] have proposed a modified Lemaitre damage model. They defined a threshold on the rate of release of the strain energy density. Khelifa et al. [28] evaluated a coupling between anisotropic elastoplastic behavior with mixed nonlinear isotropic and kinematic hardening and an isotropic ductile damage. A deep drawing test shows the efficiency of this model.

This paper discusses an improved elastoplastic model that is strongly coupled to the effect of damage in the yield function and plastic potential. The plastic flow anisotropy was enhanced by using quadratic Hill 1948 criterion based on plastic strain evolution. A mathematical formulation of the plasticity fully coupled with the damage model is developed in order to study stress changes and damage during forming operations. Moreover, the fully coupled anisotropic elastoplastic/damage model presented in this paper leads to only two nonlinear local scalar equations that are solved by the using Newton method. The constitutive equations are related to anisotropic plastic formulation, taking into account the isotropic/kinematic hardening and coupled to ductile damage model. After that, a numerical simulation of the forming process is presented. The developed model is applied in the context of the deep drawing operation using ABAQUS software via VUMAT and UMAT subroutines.

2. Plasticity Coupled with Damage

In this section, a thermodynamic model is developed to correctly describe the performances of the materials in forming processes and to develop the consistent user-friendly numerical tools. In addition, numerical models are able to capture numerous phenomena that arise through plastic deformation such as anisotropic yielding, nonlinear isotropic hardening, kinematic hardening, and damage. In this context, the elastoplastic equations defined below take into account the anisotropic plasticity and the mixed isotropic/kinematic hardening and damage.

2.1. State Variables

The simplest practice used to define the thermodynamic state of the material is the use of measured variables capable at a given point to describe the state of the material and to expect its evolution. These variables are the effective state variables that define the fictitious undamaged state of the material. From a phenomenological point of view, the concept underlying the present model is related to the subsequent pairs of effective state variables:

-

The effective stress tensors and elastic strain $\left(\tilde{\sigma },{\tilde{\epsilon }}^e\right)$;

-

The effective isotropic hardening parameters $\left(\tilde{r},\tilde{R}\right)$;

-

The effective kinematics hardening parameters $\left(\tilde{\alpha },\tilde{\mathit{X}}\right)$.

These variables can be expressed as follows

$${\tilde{\epsilon }}^e=\sqrt{1-d}{\epsilon }^e,\tilde{\sigma }=\sigma /\sqrt{1-d}$$

(1)

$$\tilde{r}=\sqrt{1-d}r,\tilde{R}=R/\sqrt{1-d}$$

(2)

$$\tilde{\alpha }=\sqrt{1-d}\alpha ,\tilde{\mathit{X}}=\mathit{X}/\sqrt{1-d}$$

(3)

where d is the isotropic ductile damage variable, R is the drag stress in isotropic hardening, $\alpha$ and r are internal variables corresponding to kinematic and isotropic hardening, respectively.

2.2. Fundamental Equations

In the incremental plasticity theory, the strain tensor can be expressed in the additive form of the elastic-damage ${\epsilon }^e$ and the plastic-damage ${\epsilon }^p$

$${\epsilon }^{}={\epsilon }^e+{\epsilon }^p$$

(4)

The Helmholtz free energy is obtained using the following form

$$\psi \left({\epsilon }^e,{\alpha }_k,r,d\right)={\psi }^e\left({\epsilon }^e,d\right)+{\psi }^p\left({\alpha }_k,r,d\right)$$

(5)

In which the elastic and plastic parts are defined by

$$\begin{gathered} {\psi }^e\left({\epsilon }^e,d\right)=\frac{1-d}{2}{\epsilon }^e:\mathit{D}:{\epsilon }^e, \\ {\psi }^p\left({\alpha }_k,r_k,d\right)=\frac{1-d}{2}\left({\displaystyle \sum _{k=1}^N{Q}_k{r}_k^2}+{\displaystyle \sum _{k=1}^M{a}_k{\alpha }_k:{\alpha }_k}\right) \end{gathered}$$

(6)

where N and M are the number of variables corresponding to isotropic and kinematic hardening, respectively, $Q_k$ and $a_k$ are material parameters, and $\mathit{D}$ is the general elastic operator of the undamaged material. The associated thermodynamic variables can be formulated from the following relationships

$$\sigma =\frac{\partial {\psi }^e}{\partial {\epsilon }^e},{\mathit{X}}_k=\frac{\partial {\psi }^p}{\partial {\alpha }_k},R=\frac{\partial {\psi }^p}{\partial r},Y=-\frac{\partial \psi }{\partial d}$$

(7)

where Y is the associated damage variable. Using Equations (6) and (7), the variables of type stress are expressed as

$$\sigma =\left(1-d\right)\mathit{D}:{\epsilon }^e,\mathit{X}={\displaystyle \sum _{k=1}^M{\mathit{X}}_k}=\left(1-d\right){\displaystyle \sum _{k=1}^M{a}_k{\alpha }_k}$$

(8)

$$\begin{gathered} R={\displaystyle \sum _{k=1}^N{R}_k}=\left(1-d\right){\displaystyle \sum _{k=1}^N{Q}_k{r}_k}, \\ Y=\frac{1}{2}{\epsilon }^e:\mathit{D}:{\epsilon }^e+\frac{1}{2}\left({\displaystyle \sum _{k=1}^N{Q}_k{r}_k^2}+{\displaystyle \sum _{k=1}^M{a}_k{\alpha }_k:{\alpha }_k}\right) \end{gathered}$$

(9)

2.3. Hill Anisotropy Criterion

The Hill 1948 model is specially suggested for anisotropic sheet metal produced using the rolling process [29]. Numerous other criteria exist in the literature [3,30] to describe the anisotropic performance of sheet metal using either quadratic or non-quadratic yield functions. The Hill criterion is the most useful in the literature due to its reduced number of parameters to be determined experimentally. The quadratic Hill 1948 criterion expresses the equivalent stress with an anisotropy operator of six parameters ($F,G,H,N,M$ and $L$) as

$${\sigma }_{eq}^H(\xi )={\left[H{\left({\xi }_1-{\xi }_2\right)}^2+F{\left({\xi }_2-{\xi }_3\right)}^2+G{\left({\xi }_3-{\xi }_1\right)}^2+2N{\xi }_{12}^2+2L{\xi }_{23}^2+2M{\xi }_{31}^2\right]}^{1/2}$$

(10)

where $F,G,H,N,M$ and $L$ are material constants obtained by tests of the material in different orientations and $\xi$ is the effective stress tensor, given by

$${\xi }_{}=\sigma -\mathit{X}$$

(11)

The quadratic criterion of Hill 1948 is integrated into the model using the fourth rank tensor $\mathit{P}$ defined by:

$${\sigma }_{eq}^H\left(\xi \right)={\phi }_f\left(\xi \right)={\Vert \xi \Vert }_{\mathit{P}}=\sqrt{{\xi }^t\mathit{P}\xi }$$

(12)

$$\mathit{P}=\left[\begin{array}{cccccc}H+G& -H& -G& 0& 0& 0\\ & H+F& -F& 0& 0& 0\\ & & F+G& 0& 0& 0\\ & & & 2N& 0& 0\\ & Sym& & & 2M& 0\\ & & & & & 2L\end{array}\right]$$

(13)

The yield and the plastic potential functions, $f$ and $F$, are given in terms of the associated thermodynamic variables, as follows

$$f=\frac{{\phi }_f\left(\xi \right)-R}{\sqrt{1-d}}-{\sigma }_Y$$

(14)

$$F=f+\frac{1}{2\left(1-d\right)}\left[{\displaystyle \sum _{k=1}^N\frac{{\beta }_k}{Q_k}{R}_k^2}+{\displaystyle \sum _{k=1}^M\frac{b_k}{a_k}{\mathit{X}}_k:{\mathit{X}}_k}\right]+F_d\left(Y,d\right)$$

(15)

where ${\sigma }_Y$ is the initial yield stress, ${\beta }_k$ and $b_k$ denote the isotropic and kinematic hardening material parameters, respectively, and $F_d\left(Y,d\right)$ can be any damage model. The evolution equations can be obtained using the generalized normal rule as follows

$${\dot{\epsilon }}^p=\dot{\gamma }\frac{\partial F}{\partial \sigma }=\frac{\dot{\gamma }}{\sqrt{1-d}}\mathit{n},\mathit{n}=\frac{1}{{\phi }_f}\mathit{P}\xi$$

(16)

$${\dot{\alpha }}_k=-\dot{\gamma }\frac{\partial F}{\partial {\mathit{X}}_k}={\dot{\epsilon }}^p-\frac{\dot{\gamma }}{\left(1-d\right)}\frac{b_k}{a_k}{\mathit{X}}_k={\dot{\epsilon }}^p-\dot{\gamma }{b}_k{\alpha }_k$$

(17)

$${\dot{r}}_k=-\dot{\gamma }\frac{\partial F}{\partial R_k}=\dot{\gamma }\left[\frac{1}{\sqrt{1-d}}-\frac{{\beta }_k{R}_k}{Q_k\left(1-d\right)}\right]=\dot{\gamma }{h}_k,h_k=\frac{1}{\sqrt{1-d}}-{\beta }_k{r}_k$$

(18)

$$\dot{d}=\dot{\gamma }\frac{\partial F}{\partial Y}=\dot{\gamma }\overline{Y},\overline{Y}=\frac{\partial F_d}{\partial Y}$$

(19)

where $\dot{\gamma }$ designates the plastic multiplier which is consistent with the Kuhn–Tucker conditions for loading and unloading

$$\dot{\gamma }\ge 0,f\le 0,\dot{\gamma }f=0$$

(20)

Using Equations (8) and (17)–(19), the evolution equation of the back stress ${\mathit{X}}_k$ is given as below:

$${\dot{\mathit{X}}}_k=\left(1-d\right)a_k{\dot{\alpha }}_k-a_k{\alpha }_k\dot{d}=\left(1-d\right)a_k\left[{\dot{\epsilon }}^p-\frac{\dot{\gamma }}{\left(1-d\right)}\frac{b_k}{a_k}{\mathit{X}}_k\right]-\frac{\dot{d}}{1-d}{\mathit{X}}_k$$

(21)

$${\dot{\mathit{X}}}_k=a_k\dot{\gamma }\sqrt{1-d}\mathit{n}-\dot{\gamma }\left(b_k+\frac{\overline{Y}}{1-d}\right){\mathit{X}}_k$$

(22)

$$\dot{\mathit{X}}={\displaystyle \sum _{k=1}^M{\dot{\mathit{X}}}_k}=\dot{\gamma }{\mathit{H}}_X,{\mathit{H}}_X=a\sqrt{1-d}\mathit{n}-\frac{\overline{Y}}{1-d}\mathit{X}-{\displaystyle \sum _{k=1}^M{b}_k{\mathit{X}}_k},a={\displaystyle \sum _{k=1}^M{a}_k}$$

(23)

3. Numerical Integration Schema

An algorithm is developed to simulate the present plasticity model. It is based on the fully implicit backward Euler integration. Indeed, for elastoplastic numerical problems, this integration scheme is usually used because of its unconditional stability. With this numerical integration procedure, it is crucial to develop the algorithmic tangent module. The main objective of using this module is to preserve the quadratic rate of the asymptotic convergence of the Newton method.

It can be followed that the enforcement of the consistency condition and the damage equation is reduced to two scalar equations. The unknowns of this system of equations are the plastic multiplier $\Delta \gamma$ and the damage variable $d_{n+1}=d$. The Newton–Raphson method is used to solve this system of equations. Appendix A summarizes the one-step to answer the couple $\left(\Delta \gamma ,d_{n+1}=d\right)$.

The derivation of the consistent tangent modulus follows only the standard application of consistent linearization concepts. In Appendix B, we illustrate its final expression. The tangents modulus for the present coupled anisotropic plasticity-ductile damage is appropriate for plane strain, plane stress, and axisymmetric and 3D problems. Only two terms are modified, which are D and P.

Finally, the complete and fully implicit integration scheme of the coupled plastic-damage with nonlinear isotropic/kinematic hardening is shown in

Table 1. Integration scheme of the coupled plastic-damage (UMAT).

Compute trial elastic stress (Elastic predictor)      $\Delta \gamma =0$      ${\sigma }^{trial}={\sigma }_n+\left(1-d_n\right)\mathit{D}\Delta \epsilon ,$ ${\overline{\mathit{X}}}_{k,n}={\mathit{X}}_{k,n}$ $\overline{\xi }={\sigma }^{trial}-{\displaystyle \sum _{k=1}^M{\overline{\mathit{X}}}_{k,n}}$ (i) Check plastic condition IF $\frac{1}{{\left(1-d\right)}^{0.5}}\left[{\phi }_{n+1}-R\right]-{\sigma }_Y\le 0$ THEN (Elastic) Set ${\left(\right)}_{n+1}={\left(\right)}_n$, ${\mathit{D}}^{ep}=\mathit{D}$ and RETURN ELSE (Plastic correction) (ii) Find $\left(\Delta \gamma ,d\right)$ by local iteration (Appendix A) (iii) Update variables      ${\xi }_{n+1}={\mathit{I}}_c^{-1}.\overline{\xi },$ ${\mathit{n}}_{n+1}=\frac{1}{{\phi }_{n+1}}\mathit{H}{\xi }_{n+1}$ ${\mathit{X}}_{k,n+1}=\frac{1-d_{n+1}}{1-d_n}{\overline{\mathit{X}}}_{k,n}+{\omega }_k{\mathit{n}}_{n+1},$ ${\mathit{X}}_{n+1}={\displaystyle \sum _{k=1}^M{\mathit{X}}_{k,n+1}},$ ${\sigma }_{n+1}={\xi }_{n+1}+{\mathit{X}}_{n+1}$ (iv) Consistent elastoplastic modulus (Appendix B)      ${\mathit{D}}^{ep}=\left(1-d_{n+1}\right)\left[{\mathit{D}}^{*}-\frac{1}{q}{\mathit{V}}_1\otimes {\mathit{V}}_2\right]$ ENDIF

.

Damage Model

For the coupled anisotropic plasticity damage models, the damage potential, $F_d$ should be considered in the derivation of the diverse equations. In this present work, the modified Lemaitre damage model is adopted and written as follows

$$F_d=\frac{S}{\left(s+1\right){\left(1-d\right)}^{\beta }}{\langle \frac{Y-Y_0}{S}\rangle }^{s+1},\overline{Y}=\frac{1}{{\left(1-d\right)}^{\beta }}{\langle \frac{Y-Y_0}{S}\rangle }^s$$

(24)

where $\beta$, $s$, $S$ and $Y_0$ are material parameters. The strain density release rate Y, defined in Equation (14), can be derived using

$$Y=Y\left({\tilde{\sigma }}_{n+1},r_k,{\alpha }_k\right)=\frac{{\tilde{q}}^2}{6G}+\frac{{\tilde{p}}^2}{2K}+\frac{1}{2}\left({\displaystyle \sum _{k=1}^N{Q}_k{r}_k^2}+{\displaystyle \sum _{k=1}^M{a}_k{\alpha }_k:{\alpha }_k}\right)$$

(25)

where K and G are the bulk and shear modulus, respectively, and $\tilde{p}$ and $\tilde{q}$ are expressed as

$$\tilde{p}=\frac{1}{3}tr\left({\tilde{\sigma }}_{n+1}\right),\tilde{q}=\sqrt{3J_2\left({\tilde{\sigma }}_{n+1}\right)}$$

(26)

4. Numerical Results of Coupled Anisotropic Plasticity and Damage Models

The main objective of this section is to develop a numerical analysis in order to validate the coupled anisotropic plasticity and damage models. The isotropic hardening describes the sheet material’s behavior during a cold forming operation. It determines the size of the yield surface. The isotropic hardening law depends on the initial value of yield stress ${\sigma }_y$ and material constants Q and β, which represent the isotropic hardening parameters. This law is given in Equation (27).

$$\sigma ={\sigma }_y+Q(1-e^{-\beta .{\epsilon }^p})$$

(27)

Previous sections describe the developed numerical models, which are implemented in the finite element software Abaqus/Explicit through UMAT and VUMAT subroutines. Four numerical examples are analyzed in the next parts.

4.1. Perforated Square Plate under Biaxial Extension

As shown in Figure 1a, the square plate has as a thickness of t = 1.0 mm and a central circular hole. The damage field is illustrated with different values of anisotropic parameter α.

The material properties are given in

Table 2. Material properties of perforated square steel plate.

Elasticity	$E=210\mathrm{GPa}$, $\nu =0.3$
Isotropic Hardening (MPa)	$Q=3300\mathrm{MPa}$; ${\sigma }_y=620\mathrm{MPa}$; $\beta =0.4$
Damage	$\delta =1$, $\beta =1$, $s=1$, $S=3.5\mathrm{MPa}$, $Y_0=0\mathrm{MPa}$

. Only isotropic hardening is considered in this example. $\alpha$ is a parameter that controls the anisotropy.

Anisotropic Hill criterion is considered with the following material properties.

$$R_{11}=1,R_{22}=\alpha R_{11},R_{33}=R_{12}=R_{13}=R_{23}=1$$

(28)

Subsequently, the classical Hill 1948 coefficients are obtained using Equation (29).

$$\left[\begin{array}{c}F\\ G\\ H\end{array}\right]=\frac{1}{2}\left[\begin{array}{ccc}-1& 1& 1\\ 1& -1& 1\\ 1& 1& -1\end{array}\right]\left[\begin{array}{c}1/R_{11}^2\\ 1/R_{22}^2\\ 1/R_{33}^2\end{array}\right],N=M=L=1.5$$

(29)

By noting the symmetry, only a quarter part of the plate is considered in the simulation. In this example, a biaxial extension is considered and applied on the top and right boundaries of the quarter plate. This plate is discretized using 4-node quadrilateral plane stress element. The finite element mesh employed is shown in Figure 1. Appropriate symmetry boundary conditions are applied. Loading is performed by controlling the same displacement on the top and right boundaries of the square plate, $u_1=u_2=u$. Depending on the choice of anisotropic parameter α, different values of the prescribed displacement $u$ are to be interpreted as ensuring that the damage does not reach an excessive value. The damage at point A (d_A) and maximum damage in the plate (d_max) are listed in

Table 3. Specified displacement and damage for different value of α.

$\mathit{\alpha } = {\mathit{R}}_{22}/{\mathit{R}}_{11}$	u (mm)	dA	dmax
0.5	0.115	0.00278	0.667
1	1.415	0.397	0.397
1.5	0.38	0.464	0.464

for different values of the parameter α.

Damage contour plots at the final state for each value of parameter α are illustrated in Figure 1. The isotropic J₂ elastoplasticity corresponds to $\alpha =1.0$. The damage contours obtained using numerical analysis are in agreement with those expected. As the parameters α increase, the damage at point A increase. For the three values of the parameter α, damage evolution at point A (d_A) versus displacement $u$ are illustrated in Figure 2.

4.2. Numerical Simulations of Bulge Test

As a second numerical test, a sheet metal forming process is analyzed in this part. In fact, the bulge test is considered. It is usually used to characterize the formability of thin sheets. This process is schematized in Figure 3.

During the bulge test, the blank lies on a rigid die and its edged is clamped. Only one quarter is modeled due to symmetry. The fully coupled anisotropic plasticity/damage model is implemented using ABAQUS/Explicit software by using the user interface subroutine VUMAT. The sheet has a thickness of 1 mm. It is meshed using 588 thin shell elements with a reduced integration of type S4R. The matrix is considered as a rigid body. It is discretized with the rigid elements of type R3D4. A hydraulic pressure is applied on the bottom surface of the circular blank sheet. The maximum value of the ramped pressure evolution is 7.2 MPa. The material properties are listed in

Table 4. Material properties considered in bulge test.

Elasticity	$E=210\mathrm{GPa}$, $\nu =0.3$
Isotropic Hardening (MPa)	${\sigma }_p=488.35{\left(0.001+{\epsilon }^p\right)}^{0.2427}$
Initial yield stress, ${\sigma }_{y_0}$(MPa)	175
Damage	$\delta =1$, $\beta =1$, $s=1$, $S=1.25\mathrm{MPa}$, $Y_0=0\mathrm{MPa}$

, according to [31]. Simulations were carried out considering an isotropic material model (r₀ = r₄₅ = r₉₀ = 1).

The damage distribution in the studied bulge test is presented in Figure 4. Based on the numerical results (Figure 4), the isotropic model shows that the maximum value of damage (SDV4) is located in the bulge center. The damaged element deletion method is activated in the simulation of bulge test using ABAQUS/Explicit. Figure 4 shows that, if the ductile damage reaches a critical value (0.98), then elements are supposed fully damaged and detached from the structure. In the same context, Figure 5 illustrates the damage evolution for the isotropic case.

The numerical damage evolutions are examined by using the nodes situated in the middle layer. Therefore, and as shown is Figure 5, the result obtained by the present model was compared with that developed by [31]. A good agreement between both results is showed.

4.3. Model Evaluation Using Uniaxial Tensile Test

In this section, the predictive model capability is examined by comparing the numerical and experimental results obtained from the uniaxial tensile test. The anisotropic elastoplastic mechanical model of DD13 sheet steel is considered in the simulation of the tensile test. The DD13 steel material is hot rolled sheet metal. This sheet metal has a ferritic microstructure. The chemical composition of the blank is presented in

Table 5. Chemical composition of the DD13 sheet material.

C	Smax	Mnmax	Pmax
0.08	0.03	0.4	0.03

. The material properties of the blank sheet are provided in

Table 6. Mechanical properties of the DD13 material [32].

Elastic Prop.			Hill 1948 Coefficients
$E\left(\mathrm{GPa}\right)$	$\nu$	${\sigma }_Y(\mathrm{MPa})$	F	G	H	N
220	0.3	278	0.43	0.47	0.53	1.49

and

Table 7. Hardening and damage parameters for DD13 [32].

Isotropic hardening parameters	$Q_1$ (MPa)	$Q_2$ (MPa)	${\beta }_1$	${\beta }_2$
	490	50	1.35	30
Damage parameters	$\beta$	$s$	$S$ (MPa)	$Y_0$ (MPa)
	5	1.5	200	15

. A double exponential law is adopted for the anisotropic elastoplastic mechanical model. The expression of the stress evolution is written as

$$\sigma \left(p\right)={\sigma }_Y+A_1\left(1-e^{{\beta }_1{\epsilon }^p}\right)+A_2\left(1-e^{{\beta }_2{\epsilon }^p}\right)$$

(30)

where ${\epsilon }_p$ represents the equivalent plastic strain, ${\sigma }_Y$ is the yield stress, and $A_k=\frac{Q_k}{{\beta }_k}$, ${\beta }_k,k=1,2$ are the material parameters.

The eight nodes hexahedral elements (C3D8) mesh is adopted for the simulation of the tensile test, the results of which is illustrated in Figure 6a. An imposed displacement was applied in the axial direction. Numerical and experimental results of the force–displacement curves of the tensile test along the rolling direction (RD) are obtained using the VUMAT subroutine and are illustrated in Figure 7.

Figure 6b presents the specimen damage for 21 mm displacement in the RD using the present model. The experimental specimen final fracture in the RD is given in Figure 6c. As depicted experimentally, the numerical damage takes place at the central zone, and the rupture facies are not perpendicular to the loading axis. Results in terms of the load–displacement curves are presented in Figure 7. By using the presented anisotropic plasticity/damage model, Figure 7 shows an acceptable fit of the numerical result when compared to experimental one. Next, the efficiency of the presented model will be evaluated in the case of the cross-die deep drawing test.

4.4. Numerical Simulation of Cross-Die Deep Drawing Test

Cross-die deep drawing simulation is carried out in this section. This test is used to cover a wide range of stress states and to reproduce most significant strain paths that can be found in deep drawing, i.e., tensile, shear, plane strain, and biaxial loading modes. An improved elastoplastic model strongly coupled with damage effect in the yield function and plastic potential is used. The developed model is implemented in ABAQUS software via the VUMAT subroutine. For the simulation, a 3D finite element model is performed in Abaqus software to analyze the cross-die test. The die, punch, and blank holder are defined as discrete rigid parts, and the sheet was assigned as a deformable homogenous shell. The initial dimensions of the square blank sheet were equal to 260 × 260 × 1 mm³. A four nodes quadrilateral shell element S4R with a size of 5 × 5 mm² and with five integration points through the thickness was employed for meshing the sheet. Penalty contact is used in forming simulations. The Coulomb friction coefficient of 0.05 was considered between sheet and tools surfaces. The blank holding force during forming was 450 kN. The punch stroke of 60 mm was considered in the simulation. Figure 8 shows the main geometrical parameters of the tools. The die and punch radii are 20 mm and 14 mm, respectively. In all cases, the blank sheet orientation with respect to the tools is shown in Figure 8. In the simulation of cross-die deep drawing, the blank is made of DD13 sheet metal.

After the simulation, the blank sheet thinning was determined along three directions (i.e., rolling direction (RD), transverse direction (TD), and diagonal direction (DD)) as illustrated in Figure 9. From this figure, it is noticed that the thinning of the simulated workpiece is more important in DD (0.49 mm-thickness) compared to RD and TD (≈0.7 mm-thickness).

Figure 10 and Figure 11 show the simulation results of the cross-die test. The damage location and the force–displacement response were predicted precisely, as presented in the work of [33,34].

The minor and major strains of the 60 mm-height of the cross-die are illustrated in Figure 12. The strain path analysis of the cross-die simulation reveals the complex forming process. In the simulated part, the strain distributions of some points are near or above the experimental forming limit diagram (FLD) curve. The location of the damage is the same as that presented in the experimental work of [33]. In other words, the figure can illustrate the capability of the proposed model in the damage prediction of a cross-die deep drawing test.

5. Conclusions

In this paper, a coupled model of anisotropic plasticity and continuous ductile damage is proposed. A generalized quadratic yield function is considered to exhibit isotropic and nonlinear kinematic hardening. According to the literature, more than two local equations are required in coupled elastoplastic damage models. However, in our work, only two scalar equations are required to solve the developed numerical algorithm. The unknowns in this system of equations are the plastic multiplier and the damage variables. Using the implicit solver, and in order to preserve the quadratic rate of asymptotic convergence of Newton’s method, the algorithmic tangent modulus is developed. It is given in closed form using exact linearization. The mathematical models are implemented using the ABAQUS/ Explicit software and the UMAT and VUMAT subroutines of the user interface. The results of DD13 sheet metal in cross-die deep drawing show the effect of anisotropic coefficients and hardening and damage parameters on the distribution of thickness and damage in the obtained workpiece. Numerical examples are developed which show the efficiency and ability of the proposed model to predict ductile damage growth in the sheet forming process.