Prediction of Plastic Instability in Sheet Metals During Forming Processes Using the Loss of Ellipticity Approach

Abstract

The prediction of plastic instability in sheet metals during forming processes represents nowadays an ambitious challenge. To reach this goal, a new numerical approach, based on the loss of ellipticity criterion, is proposed in the present contribution. A polycrystalline model is implemented as a user-material subroutine into the ABAQUS/Implicit finite element (FE) code. The polycrystalline constitutive model is assigned to each integration point of the FE mesh. To derive the mechanical behavior of this polycrystalline aggregate from the behavior of its microscopic constituents, the multiscale self-consistent scheme is used. The mechanical behavior of the single crystals is described by a finite strain rate-independent constitutive framework, where the Schmid law is used to model the plastic flow. The condition of loss of ellipticity at the macroscale is used as plastic instability criterion in the FE modeling. This numerical approach, which couples the FE method with the self-consistent scheme, is used to simulate a deep drawing process, and the above criterion is used to predict the formability limit of the studied sheets during this operation.

Keywords:
crystal plasticity; loss of ellipticity; self-consistent; finite elements; forming limit diagrams

1 INTRODUCTION

For weight reduction of automotive parts, there is nowadays a growing need for designing new lightweight high strength metallic materials (Heller et al., 1998). To meet this industrial demand, these new materials should present some desirable in-use properties, such as high strength and improved ductility. They are generally characterized by their strong anisotropy, which is the most important physical phenomenon that affects the strength and the formability of sheet metals (Guan et al., 2006). Considering the high costs relating to the experimental characterization of these new alloys, numerical simulation tools are expected to provide a valuable alternative to be used to assess the strength and the formability of these new materials. This performance evaluation is considered to be an essential step in the strategy of the design of industrial parts produced from these materials (Nakamachi, 1997). In order to evaluate their formability during sheet metal forming processes, finite element (FE) analyses are frequently used in the automotive industry. However, reliable FE simulations of forming processes require an accurate description of the initial and induced anisotropy of the studied sheet metals. Unfortunately, phenomenological approaches, which have been widely used for several decades in the FE simulations, are not able to accurately account for this initial and induced anisotropy. To overcome this limitation, various micromechanical models have been developed and incorporated into FE analyses to simulate different forming processes (Beaudoin et al., 1994; Guan et al., 2006; Nakamachi et al., 2001). The advantage of such micromechanical modeling, compared to phenomenological approaches, is its ability to link, in a natural way, several material features (initial and induced anisotropy, grain morphology…) to some important in-use properties (strength, formability …). Thanks to the substantial increase in computational power, the use of micromechanical constitutive frameworks in the simulation of complex forming processes becomes an interesting and feasible alternative to phenomenological models, the latter having been intensively used in the past.

In the present investigation, the numerical strategy based on the coupling between micromechanical modeling and FE analysis is followed. The self-consistent approach (Lipinski and Berveiller, 1989) is used as multiscale scheme in this micromechanical modeling. Within this approach, each single crystal is assumed to be an ellipsoidal inclusion contained in a fictitious homogeneous medium endowed with the same mechanical features as the polycrystalline aggregate. In contrast to the Taylor model, which is more commonly used in this kind of coupled strategies (Guan et al., 2006; Jung et al., 2013), the equilibrium equations at the grain level are fully satisfied when the self-consistent approach is used. Moreover, this choice allows considering some effects that cannot be taken into account with the Taylor model, such as the impact of grain morphology. Furthermore, the Taylor model may be viewed as a particular case of the self-consistent approach. The mechanical behavior of the single crystals is described by a finite strain rate-independent constitutive framework, where the Schmid law is used to model the plastic flow.

The associated self-consistent approach has been fully implemented, as a user-material (UMAT) subroutine, into the ABAQUS/Implicit FE code. To predict the occurrence of plastic instability in the studied sheet metals during forming processes, a theoretical criterion, based on the loss of ellipticity of the macroscopic governing equilibrium equations (Rudnicki and Rice, 1975), has been integrated into the UMAT subroutine. Note that, in the context of localized necking in sheet metals, the application of this plastic instability criterion in conjunction with the developed rate-independent self-consistent approach allows predicting forming limit diagrams (FLDs) at realistic strain levels for all strain paths, contrary to the case of classical phenomenological models, which require introducing additional destabilising phenomena, such as vertex-like yield surface effects or softening mechanisms. The main reason to this is that typical destabilizing effects are a natural outcome of crystal plasticity, through the formation of vertices at the current points of the Schmid yield surface. The implemented UMAT subroutine has been first used within a single eight-node hexahedral element (C3D8R), with only one integration point, in order to reproduce various homogeneous stress and strain states until the condition of loss of ellipticity is satisfied. The use of a single finite element in this preliminary application allows assimilating the studied medium to an infinite block of material, which is homogeneously deformed until the incipience of localization bifurcation. Accordingly, the loss of ellipticity phenomenon can be viewed as the result of an intrinsic material instability with no structural effects. In this first application, the material instability of the studied polycrystalline aggregate has been represented through the well-known concept of forming limit diagrams (Keeler and Backofen, 1963). Then, the above-described numerical implementation has been used to simulate the deep drawing of various sheet specimens. It is noteworthy that this second application allows the evolution of microstructural properties (crystallographic texture as well as grain morphology) to be accurately predicted during deep drawing processes. In these deep drawing simulations, the loss of ellipticity phenomenon may be attributed to competing material and structural instabilities.

The reminding parts of the paper are structured as follows:

• The second Section is devoted to the presentation of the constitutive frameworks both at the single crystal and polycrystal scales.

• The theoretical framework for the loss of ellipticity approach is presented in the third Section.

• In the fourth Section, two numerical applications mainly are shown to illustrate the efficiency of the implemented numerical strategy.

Notations, conventions and abbreviations

The notations, conventions and abbreviations used in this paper are clarified in the box bellow. Additional notations will be provided as needed.

Vector and tensor fields are designated by bold letters and symbols. Scalar variables and parameters are represented by thin letters and symbols. Microscopic (resp. macroscopic) variables are referenced by lowercase (resp. capital) letters. Einstein’s convention of implied summation over repeated indices is adopted. The range of free (resp. dummy) index is given before (resp. after) the corresponding equation. • ˙ time derivative of • I 2 the second-order identity tensor I 4 the fourth-order identity tensor •∇ lattice co-rotational derivative of • •-1 inverse of • •. • inner product • : • double contraction product • ⊗ • tensor product •T transpose of second-order tensor • det(•) determinant of second-order tensor • Tr(•) trace of second-order tensor • sgn(•) sign of •

Vector and tensor fields are designated by bold letters and symbols.
Scalar variables and parameters are represented by thin letters and symbols.
Microscopic (resp. macroscopic) variables are referenced by lowercase (resp. capital) letters.
Einstein’s convention of implied summation over repeated indices is adopted. The range of free (resp. dummy) index is given before (resp. after) the corresponding equation.
$$\dot{•}$$	time derivative of •
I 2	the second-order identity tensor
I 4	the fourth-order identity tensor
•∇	lattice co-rotational derivative of •
•-1	inverse of •
•. •	inner product
• : •	double contraction product
• ⊗ •	tensor product
•T	transpose of second-order tensor •
det(•)	determinant of second-order tensor •
Tr(•)	trace of second-order tensor •
sgn(•)	sign of •

2 CONSTITUTIVE FRAMEWORKS AT THE SINGLE CRYSTAL AND POLYCRYSTAL SCALES

In this section, several developments will be carried out in order to express the behavior laws both at the single crystal and polycrystal scales in the following generic forms:

$$\begin{array}{l}\text{single crystal behavior law (microscopic level):}\dot{n}=l:g\\ \text{polycrystal behavior law (macroscopic level):}\dot{N}=\mathcal{L}:G,\end{array}$$ (1)

where:

• $$\dot{n}$$ (resp. $$\dot{N}$$) is the microscopic (resp. macroscopic) nominal stress rate tensor.

• l (resp. ℒ) is the microscopic (resp. macroscopic) tangent modulus.

• g (resp. G) is the microscopic (resp. macroscopic) velocity gradient tensor.

The aim of the following sections is to derive the explicit expressions for both the microscopic (Section 2.1) and macroscopic (Section 2.2) tangent moduli.

2.1 Single Crystal Scale

The microscopic velocity gradient g can be additively decomposed into its symmetric and skew-symmetric parts denoted as d and w, respectively

$$g=d+w$$ (2)

Considering the fact that elastic deformation is sufficiently small compared to plastic deformation, d and w can be in turn split into their elastic and plastic parts

$$d=d^e+d^p;w=w^e+w^p$$ (3)

The rotation r of the single crystal lattice frame is determined by we and satisfies the relation below

$$\dot{r}.r^T=w^e$$ (4)

Because plastic slip occurs along specific crystallographic directions, on particular crystallographic planes, the components _^dp and _^wp of the plastic part _^gp of the velocity gradient can be written in the following form:

$$\{\begin{array}{c}{d}^p={\dot{\gamma }}^{\alpha }\mathrm{sgn}({\tau }^{\alpha })R^{\alpha }\\ w^p={\dot{\gamma }}^{\alpha }\mathrm{sgn}({\tau }^{\alpha })P^{\alpha }\end{array};\alpha =\mathrm{1,...,}{N}_s$$ (5)

where:

$${\dot{\gamma }}^{\alpha }$$ is the absolute value of the slip rate on the (-th slip system.

_^Ra and _^Pa are the symmetric and skew-symmetric parts of the Schmid orientation tensor $${\overrightarrow{m}}^{\alpha }\otimes {\overrightarrow{n}}^{\alpha }$$ (Schmid and Boas, 1935), respectively.

$${\overrightarrow{m}}^{\alpha }$$ and $${\overrightarrow{n}}^{\alpha }$$ are the slip direction and the slip plane normal vector, respectively. In the single crystal lattice frame, the counterpart $${\overrightarrow{m}}_0^{\alpha }$$ (resp. $${\overrightarrow{n}}_0^{\alpha }$$) of $${\overrightarrow{m}}^{\alpha }$$ (resp. $${\overrightarrow{n}}^{\alpha }$$) remains fixed during the deformation. Vectors $${\overrightarrow{m}}_0^{\alpha }$$ and $${\overrightarrow{n}}_0^{\alpha }$$ are related to $${\overrightarrow{m}}^{\alpha }$$ and $${\overrightarrow{n}}^{\alpha }$$ by the following relations:

$${\overrightarrow{m}}^{\alpha }=r.{\overrightarrow{m}}_0^{\alpha };{\overrightarrow{n}}^{\alpha }={\overrightarrow{n}}_0^{\alpha }.r^T$$ (6)

_^Ns is the total number of slip systems.

τa is the resolved shear stress defined as the projection of the Cauchy stress tensor s on the Schmid orientation tensor $${\overrightarrow{m}}^{\alpha }\otimes {\overrightarrow{n}}^{\alpha }$$:

$${\tau }^{\alpha }=\sigma :({\overrightarrow{m}}^{\alpha }\otimes {\overrightarrow{n}}^{\alpha })$$ (7)

The elastic part of the behavior law classically relates the co-rotational derivative _^s∆ of the Kirchhoff stress tensor to the elastic strain rate _^de as follows:

$$s^{\nabla }={\mathcal{C}}^e:d^e$$ (8)

where _^Ce is the fourth-order elasticity tensor. Here, linear isotropic elasticity is assumed. We recall that the Kirchhoff stress tensor s is related to the Cauchy stress tensor σ by the following relation:

$$s=j\sigma$$ (9)

where j is the jacobian of the deformation gradient.

The relation between the time derivative $$\dot{s}$$ of the Kirchhoff stress tensor and its co-rotational derivative reads

$$\dot{s}=s^{\nabla }-s.w^e+w^e.s$$ (10)

The nominal stress tensor n is defined by the relation n = j f ^-1 .σ and, adopting an updated Lagrangian approach, its rate $$\dot{n}$$ reads

$$\dot{n}=\dot{\sigma }+\sigma Tr(d)-g.\sigma$$ (11)

The latter can be related to the slip rates $${\dot{\gamma }}^{\alpha }$$ by combining Eqs. (3), (5), (8)-(11)

$$\dot{n}={\mathcal{C}}^e:d-\sigma .w-d.\sigma -({\mathcal{C}}^e:R^{\alpha }+P^{\alpha }.\sigma -\sigma .P^{\alpha }){\dot{\gamma }}^{\alpha }\mathrm{sgn}({\tau }^{\alpha })$$ (12)

The plastic flow rule used to determine the expression of the slip rates $${\dot{\gamma }}^{\alpha }$$ results from the Schmid criterion (Schmid and Boas, 1935). This criterion states that plastic slip occurs only when the absolute value of the resolved shear stress _^τα on slip system α reaches a critical value $${\tau }_c^{\alpha }$$, which can be summarized in the following form:

$$\forall \alpha =\mathrm{1,...,}{N}_s:\left(\left|{\tau }^{\alpha }\right|-{\tau }_c^{\alpha }\right){\dot{\gamma }}^{\alpha }=0$$ (13)

The overall strain hardening characteristic of the material is modeled by the introduction of a hardening matrix h, which describes the evolution of the critical resolved shear stress $${\tau }_c^{\alpha }$$.

$$\forall \alpha =\mathrm{1,...,}{N}_s:{\dot{\tau }}_c^{\alpha }=h^{\alpha \beta }{\dot{\gamma }}^{\beta };\beta =\mathrm{1,...,}{N}_s$$ (14)

The following consistency condition can be derived from the Schmid law (13)

$$\forall \alpha \in \mathcal{A}:{\dot{\tau }}^{\alpha }\mathrm{sgn}({\tau }^{\alpha })-{\dot{\tau }}_c^{\alpha }=0;{\dot{\gamma }}^{\alpha }>0$$ (15)

where 𝒜 is the set of active slip systems. The slip rates of the slip systems that do not belong to set 𝒜 are obviously equal to zero.

The time derivative of the resolved shear stress $${\dot{\tau }}^{\alpha }$$ is given by the following expression:

$${\dot{\tau }}^{\alpha }={\sigma }^{\nabla }:({\overrightarrow{m}}^{\alpha }\otimes {\overrightarrow{n}}^{\alpha })$$ (16)

where _^σ∆ is defined by the following relation:

$${\sigma }^{\nabla }=\dot{\sigma }+\sigma .w^e-w^e.\sigma$$ (17)

By using Eqs. (3) and (8), $${\dot{\tau }}^{\alpha }$$ can be further developed as follows:

$$\forall \alpha =\mathrm{1,...,}{N}_s:{\dot{\tau }}^{\alpha }=R^{\alpha }:\left[{\mathcal{C}}^e:d-\sigma Tr(d)\right]-\mathrm{sgn}({\tau }^{\beta }){\dot{\gamma }}^{\beta }{R}^{\alpha }:{\mathcal{C}}^e:R^{\beta };\beta =\mathrm{1,...,}{N}_s$$ (18)

The combination of Eqs. (14), (15) and (18) leads to the following expression for the slip rates $${\dot{\gamma }}^{\alpha }$$:

$$\forall \alpha \in \mathcal{A}:{\dot{\gamma }}^{\alpha }=J^{\alpha \beta }\mathrm{sgn}({\tau }^{\beta })R^{\beta }:\left({\mathcal{C}}^e-\sigma \otimes I_2\right):d;\beta \in \mathcal{A}$$ (19)

where J is the inverse of matrix M defined by the following index form:

$$\forall \alpha ,\beta \in \mathcal{A}:M^{\alpha \beta }=\left(h^{\alpha \beta }+\mathrm{sgn}({\tau }^{\alpha })\mathrm{sgn}({\tau }^{\beta })R^{\alpha }:{\mathcal{C}}^e:R^{\beta }\right)$$ (20)

By introducing relation (19) into Eq. (12), one finally obtains the following expression for the microscopic tangent modulus l:

$$\begin{array}{l}l={\mathcal{C}}^e-{}_{\sigma }{}^{1}l-{}_{\sigma }{}^{2}l\\ -\mathrm{sgn}({\tau }^{\alpha })\mathrm{sgn}({\tau }^{\beta })J^{\alpha \beta }\left({\mathcal{C}}^e:R^{\alpha }+P^{\alpha }.\sigma -\sigma .P^{\alpha }\right)\otimes \left(R^{\beta }:\left({\mathcal{C}}^e-\sigma \otimes I_2\right)\right);\alpha ,\beta \in \mathcal{A}\end{array}$$ (21)

where the fourth-order tensors $${}_{\sigma }{}^{1}l$$ and $${}_{\sigma }{}^{2}l$$ are the terms induced by lattice rotation within the framework of large strains. Their expressions read

$${}_{\sigma }{}^{1}l{}_{ijkl}=\frac{1}{2}({\delta }_{jl}{\sigma }_{ik}-{\delta }_{jk}{\sigma }_{il});{}_{\sigma }{}^{2}l{}_{ijkl}=\frac{1}{2}({\delta }_{ik}{\sigma }_{jl}+{\delta }_{il}{\sigma }_{jk})$$ (22)

2.2 Polycrystalline Scale

The self-consistent scale-transition scheme is adopted for deriving the overall macroscopic behavior from knowledge of the behavior of individual single crystals. The main lines of the approach are recalled hereafter; the reader may refer to (Lipinski and Berveiller, 1989; Lipinski et al., 1995) for more details. The macroscopic tensor fields are defined as the volume averages of their microscopic counterparts, and they are related through fourth-order concentration tensors A and B in order to solve the averaging problem

$$g(x)=A(x):G;\dot{n}(x)=B(x):\dot{N}\text{.}$$ (23)

Making use of Eqs. (1) and (23), the generic form of the macroscopic tangent modulus ℒ can be expressed as follows:

$$\mathcal{L}=\frac{1}{V}{\displaystyle {\int }_{V}l(x):A(x)dv}=〈l(x):A(x)〉$$ (24)

For a polycrystalline aggregate composed of ellipsoidal crystals with different crystallographic orientations, the behavior and mechanical fields of each individual grain are assumed to be homogeneous. Denoting by gI (respectively lI) the volume average of the velocity gradient (respectively tangent modulus) for a grain I within the polycrystalline aggregate, and using some elaborate derivations based on Green’s tensors (Lipinski et al., 1995), we can demonstrate that the concentration tensor AI corresponding to grain I is given by

$$A^I={\left(I_4-T^{II}:(l^I-\mathcal{L})\right)}^{-1}:{〈{\left(I_4-T^{II}:(l^I-\mathcal{L})\right)}^{-1}〉}^{-1}$$ (25)

where _^TII is the interaction tensor for grain I, related to Eshelby’s tensor (Eshelby, 1957) for an ellipsoidal inhomogeneity.

Furthermore, in order to simplify the integral equation, the assumption of homogeneous deformation inside each individual grain is made. Accordingly, if we introduce the indicator function θI(x) defined as

$$\{\begin{array}{l}{\theta }^{\text{I}}(x)=\text{1 if}x\in V^{\text{I}}\\ {\theta }^{\text{I}}(x)=\text{0 if}x\notin V^{\text{I}}\end{array}$$ (26)

then

$$g(x)=g^I{\theta }^I(x);l(x)=l^I{\theta }^I(x);I=\mathrm{1,...,}{N}_g$$ (27)

where _^VI denotes the volume of grain I and _^Ng the number of grains that make up the polycrystalline aggregate.

Under the above assumption of homogeneity of the microscopic mechanical fields, Eq. (24), which defines the macroscopic tangent modulus, can be equivalently expressed as

$$\mathcal{L}=f^I{l}^I:A^I;I=\mathrm{1,...,}{N}_g$$ (28)

where _^fI represents the volume fraction of grain I.

3 LOSS OF ELLIPTICITY APPROACH

The bifurcation theory developed in Rudnicki and Rice (1975) is known to be a relevant approach for the prediction of material instabilities associated with jumps in strain rates across a localization band.

Because field equations have to be satisfied, the kinematic condition for the strain rate jump reads

$$〚G〛=\dot{\overrightarrow{C}}\otimes \overrightarrow{\mathcal{N}}$$ (29)

where:

• ⟦X⟧ is the jump of the field X across the discontinuity band.

• $$\dot{\overrightarrow{C}}$$ is the jump vector.

• $$\overrightarrow{\mathcal{N}}$$ is the unit vector normal to the localization band.

On the other hand, the continuity of the equilibrium of forces through the band is expressed as:

$$\overrightarrow{\mathcal{N}}.〚\dot{N}〛=\overrightarrow{0}$$ (30)

The combination of Eqs. (29) and (30) with the behavior law (1) leads to the following equation:

$$\overrightarrow{\mathcal{N}}.(\mathcal{L}:\dot{\overrightarrow{C}}\otimes \overrightarrow{\mathcal{N}})=\overrightarrow{0}$$ (31)

which is equivalent to

$$(\overrightarrow{\mathcal{N}}.\mathcal{L}.\overrightarrow{\mathcal{N}}).\dot{\overrightarrow{C}}=\overrightarrow{0}$$ (32)

where $$\overrightarrow{\mathcal{N}}.\mathcal{L}.\overrightarrow{\mathcal{N}}$$ is the so-called acoustic tensor. As long as this acoustic tensor is invertible, the jump vector $$\dot{\overrightarrow{C}}$$ remains equal to zero, thus precluding any discontinuity (bifurcation) in the deformation field. However, when the acoustic tensor becomes singular, there is an infinity of jump vectors that satisfy Eq. (32), and this infinity of solutions can be seen as an indicator of effective bifurcation. From a numerical point of view, strain localization occurs when the acoustic tensor is no more invertible

$$\det (\overrightarrow{\mathcal{N}}.\mathcal{L}.\overrightarrow{\mathcal{N}})=0$$ (33)

Eq. (33) is the so-called Rice bifurcation criterion (Rudnicki and Rice, 1975), which corresponds to the loss of ellipticity of the partial differential equations governing the associated boundary value problem.

4 NUMERICAL PREDICTIONS AND DISCUSSIONS

The developed polycrystalline constitutive framework and the associated loss of ellipticity criterion have been implemented into the finite element code ABAQUS/Standard as user-material (UMAT) subroutine. Various numerical simulations based on these implementations are presented in this section. After the material data used in the different simulations are specified, the implemented model is first used to predict the formability of a sheet metal in the form of forming limit diagrams. Then, a deep drawing process is simulated using the implemented model.

4.1 Material Data

In this section, the microscopic material data used to simulate the stress-strain response of polycrystalline aggregates, on the basis of the self-consistent model, are presented. The values of the material parameters are chosen in the range of those for rolled aluminum sheets. Thus, the following choices and material parameters are adopted in the remainder of the paper:

• Single crystals with FCC crystallographic structure are used in what follows. For this crystallographic structure, the number of slip systems _^Ns is equal to 12. More details on the numbering of the corresponding slip system vectors, $${\overrightarrow{m}}_0^{\alpha }$$ and $${\overrightarrow{n}}_0^{\alpha }$$, can be found in Bettaieb et al. (2012).

• The Poisson ratio and the Young modulus are set to 0.33 and 67.7 GPa, respectively.

• The initial critical resolved shear stress τ ₀ is assumed to be the same for the different slip systems and is taken equal to 40 MPa.

• Similar to several investigations conducted in the literature on FCC materials, an isotropic hardening matrix is used to describe the evolution of the critical resolved shear stresses inside the single crystals. Accordingly, the components _^hαβ of the hardening matrix h are taken to be identical and equal to h. The hardening modulus h is given by the following power law:

$$h=h_0{\left(1+\frac{h_0\Gamma }{{\tau }_{0}n}\right)}^{n-1};\Gamma ={\displaystyle \underset{0}{\overset{t}{\int }}{\displaystyle \sum _{\alpha =0}^{N_s}\left|{\dot{\gamma }}^{\alpha }\right|dt}}$$ (34)

where n is the power-law hardening exponent and h ₀ is the initial hardening rate. The values of n and h ₀ are set equal to 0.35 and 390 MPa, respectively.

In order to assess the impact of the number of grains, which constitute the polycrystalline aggregate, on the macroscopic ductility of the material, the simulations of three polycrystalline aggregates, with 50, 500, and 1000 grains, are carried out in Section 4.2. The initial crystallographic textures of these different polycrystalline aggregates are randomly generated and are displayed in Figure 1.

Figure 1
Initial texture of the studied polycrystalline aggregates, {111} pole figure: (a) 50 grains; (b) 500 grains; (c) 1000 grains.

The different polycrystalline aggregates (with 50, 500, and 1000 grains) are simulated by subjecting them to a uniaxial tension along the rolling direction. The simulated tensile stress-strain responses associated with these three polycrystalline aggregates are shown in Figure 2. It can be observed that the uniaxial tensile stress-strain responses corresponding to the polycrystals with 500 and 1000 grains are practically the same. However, the response of the polycrystal with 50 grains is relatively different from the others. These preliminary results suggest that 500 grains is a reasonable choice to be adopted in order to generate a polycrystalline aggregate that can be considered to be representative of the studied material.

Figure 2
Effect of the number of grains on the uniaxial tensile stress-strain response.

4.2 FLD Predictions

In order to numerically predict the formability of sheet metals in the form of forming limit diagrams (FLDs), a single C3D8R solid finite element is considered. As illustrated in Figure 3, displacement-type boundary conditions, which evolve over time, are prescribed to this eight-node hexahedral element. The evolution of the displacements U and V, applied to the faces CDHG and CBFG, respectively, are defined by the following relations:

$$U(t)=e^{{\dot{\epsilon }}_{11}t};V(t)=e^{\rho {\dot{\epsilon }}_{11}t}$$ (35)

where ρ is the strain-path ratio, assumed to be constant during the deformation to guarantee the proportionality of the strain paths. This strain-path ratio ρ is varied in the range -1/2≤ ρ ≤ 1 to span the complete FLD. It is noteworthy that the plane stress state is a natural outcome of the prescribed boundary conditions. Indeed, as the face EFGH is free from any constraint, the stress components Σ₁₃, Σ₂₃, and Σ₃₃ are obviously equal to 0.

Figure 3
Illustration of the boundary conditions applied to a single finite element.

The evolution of the minimum of the determinant of the acoustic tensor, over all possible localization band orientations, as a function of the macroscopic major strain E11, is shown in Figure 4 for four representative strain paths (ρ = -0.5, ρ = 0, ρ = 0.5, and ρ = 1). As the acoustic tensor is positive definite in the deformation range that precedes the occurrence of localized necking, the unit normal vector $$\overrightarrow{\mathcal{N}}$$ corresponding to the vanishing of its determinant at the onset of strain localization is the same as that minimizing this determinant. It can be seen that, irrespective of the selected strain path and the polycrystalline aggregate, the minimum of that determinant abruptly decreases during the transition between the elastic and plastic regimes; then, this determinant continues decreasing, although more moderately, until it vanishes at the occurrence of strain localization. Several phenomena accounted for by the micromechanical model are responsible for this reduction in the determinant of the acoustic tensor. Among them, one may mention the lattice rotation, the convective stress components involved in the expression of the tangent modulus, and the slip system activity governed by the Schmid law, which represents the most important factor. Indeed, with a regular form of the Schmid law, which rounds off corners at the yield surface, the strong reduction in the determinant of the acoustic tensor cannot be observed (Akpama et al., 2016). Furthermore, it is also observed from Figure 4 that the predictions given by the polycrystalline aggregate with 50 grains exhibit much more oscillations than those obtained with the other polycrystalline aggregates, which comprise a larger number of grains. Indeed, the smoothness in the simulated curves is closely related to the continuous character of the macroscopic tangent modulus, which is itself determined by the tangent moduli of the constituent grains. However, at the grain level, the tangent modulus is not smooth due to the discrete nature of the slip system activity (i.e., the set A of active slip systems introduced in Section 2.1 changes very frequently during the deformation). Consequently, a sufficiently large number of grains is required to obtain a relatively smooth evolution for the minimum of the determinant of the acoustic tensor. Moreover, Figure 4 clearly shows that the simulation results given by the polycrystal with 500 grains are quite close to those predicted by the polycrystal with 1000 grains, especially for the strain paths ρ = 0, 0.5, and 1. These results support once again the idea that a polycrystalline aggregate with 500 grains or more provides a good representation for the mechanical behavior of the polycrystalline material.

Figure 4
Evolution of the minimum of the determinant of the acoustic tensor as a function of E11 for the different polycrystalline aggregates (50, 500, and 1000 grains): (a) ρ = -0.5; (b) ρ = 0; (c) ρ = 0.5; (d) ρ = 1.

The effect of the number of grains on the shape and location of the predicted FLDs is investigated in Figure 5. This figure clearly reveals large differences between the FLD predicted by the polycrystalline aggregate with 50 grains and the other FLDs, especially in the range of positive strain paths. This result confirms the preliminary trends observed in Figure 5. Furthermore, plastic strain localization is predicted at realistic strain levels in the range of position strain-path ratios (ρ ≥ 0), in contrast to phenomenological flow theories with associative plasticity and smooth yield surface, which generally require the introduction of some softening phenomena, such as those induced by ductile damage (Mansouri et al., 2014). The realistic level for the predicted limit strains is made possible here thanks to the yield surface vertex structure (inherent in crystal plasticity theory), which arises from the discrete nature of crystalline slip.

Figure 5
Forming limit diagrams for the different polycrystalline aggregates (50, 500, and 1000 grains).

4.3 Deep Drawing Simulations

In order to efficiently simulate the deep drawing process using the developed multiscale model, a sensitivity analysis has been carried out in order to select the optimal number of grains to be used for representing the polycrystalline aggregate associated with each macroscopic integration point of the mesh. Based on this study, it is found that an aggregate with 500 grains seems to be a pragmatic compromise in terms of good representativeness of the actual material and reasonable CPU times required to run the various FE simulations. The initial crystallographic texture used for all of the integration points in the various FE models studied in the current Section is the same as that presented in Figure 1 (b).

For predicting the macroscopic instability (at the polycrystalline aggregate scale) with the proposed model, a hemispherical punch test is performed. The description of the tool set-up is shown in Figure 6. The diameter of the punch is taken equal to 100 mm. For the friction between the punch and the blank, the associated Coulomb friction coefficient is assumed to be equal to 0.1. The plane geometries and the dimensional characteristics of the used specimens are the same as those defined in Kami et al. (2015), as displayed in Figure 7. However, the thickness is taken equal to 1 mm for all of the studied specimens. Three values for the width parameter w are used in the FE simulations: 30 mm, 90 mm, and 185 mm (the latter corresponds to a circular specimen). The specimen width w is varied in order to reproduce different strain paths in the center of the specimen: the values 30 mm, 90 mm, and 185 mm correspond to uniaxial tensile state (ρ = -0.5), plane-stress tensile state (ρ = 0), and equi-biaxial tensile state (ρ = 1), respectively. The different tools (die, punch, and blank holder) are modeled as analytical rigid surfaces, while the specimens are meshed using C3D8R linear solid finite elements. In order to capture the bending effects and hence to improve the computation accuracy, three layers of elements are used in the thickness direction. This choice seems to be a pragmatic compromise in terms of good accuracy of the results and reasonable computational times. The mesh density is increased in the contact zone between the punch and the specimen. To carry out the forming process, the punch is displaced in the vertical direction up to a maximal displacement of 60 mm.

Figure 6
Schematic description of the tools used in the simulation of the deep drawing process.

Figure 7
Dimensional characteristics of the circular and notched specimens used in the Nakajima tests (Kami et al., 2015).

Different punch force versus displacement curves are shown in Figure 8. As expected, it is observed from the curves of Figure 8 that as the width w increases, there is an increase in the load requirement. The dot reported on the different curves indicates the occurrence of macroscopic instability, as predicted by the loss of ellipticity criterion within the associated specimens. It is clearly shown that the onset of such plastic instability does not necessarily coincide with the maximum of the force punch for the various specimens.

Figure 8
Punch force-displacement curves for the different specimen geometries.

Figure 9 shows the deformed shapes for the various Nakajima specimens, with the associated von Mises equivalent stress distribution, at the onset of plastic instability. One may notice that the value of equivalent stress is relatively low in the central area of the specimens. However, this von Mises equivalent stress reaches its maximal values in narrow zones near the dome apex.

Figure 9
Distribution of the von Mises equivalent stress in the various Nakajima specimens at the onset of plastic instability: (a) Specimen with w = 30 mm; (b) Specimen with w = 90 mm; (c) Specimen with w = 185 mm.

In order to illustrate the occurrence and development of macroscopic instability, the distribution of the loss of ellipticity indicator at the onset of instability is displayed in Figure 10. This indicator takes as value 0 (resp. 1) when the minimum of the determinant of the acoustic tensor, over all possible localization band orientations, is positive (resp. negative). Red (resp. blue) color is used in Figure 10 to designate the finite elements where the loss of ellipticity indicator is equal to 1 (resp. 0). The plastic instability, which corresponds to the nullity of the loss of ellipticity indicator, is first reached in the dome apex region, where the von Mises equivalent stress is maximal (see Figure 9). For more clarity, the corresponding elements in red color are surrounded by more visible red circles in Figure 10.

Figure 10
Distribution of the loss of ellipticity indicator in the different Nakajima specimens at the onset of plastic instability: (a) Specimen with w = 30 mm; (b) Specimen with w = 90 mm; (c) Specimen with w = 185 mm.

The distribution of the loss of ellipticity indicator in the different Nakajima specimens at the end of the simulation (punch depth = 60 mm) is shown in Figure 11. It is observed from this figure that the instability zone has extended to the entire dome apex region.

Figure 11
Distribution of the loss of ellipticity indicator in the different Nakajima specimens at the end of the simulation (punch depth = 60 mm): (a) Specimen with w = 30 mm; (b) Specimen with w = 90 mm; (c) Specimen with w = 185 mm.

In order to better understand the initiation and evolution of the instability, three particular elements are selected in each specimen for the following analysis: Element A (Elem. A), which corresponds to the central element of the specimen, Element B (Elem. B), which is representative of the zone near the localization area, and Element C (Elem. C), which corresponds to the element where the instability is first triggered (see Figure 12).

Figure 12
Locations of the different analyzed elements: (a) Specimen with w = 30 mm; (b) Specimen with w = 90 mm; (c) Specimen with w = 185 mm.

The evolution of the von Mises equivalent plastic strain rate $$D_{eq}^p$$ as a function of the punch displacement for the different analyzed elements and the different specimens is shown in Figure 13. It is clearly observed from this figure that Elem. A (resp. Elem. C) has the lower (resp. higher) strain rate among the three analyzed elements. Consequently, Elem. A exhibits the slower deformation, while Elem. C displays the more rapid plastic deformation, which ultimately tends to localize. This trend is common to the three specimens. Note however that the different elements do not have the same evolution for $$D_{eq}^p$$. Indeed, for Elem. A and Elem. B, $$D_{eq}^p$$ increases before the onset of instability in the specimen, while it decreases afterwards, which typically corresponds to elastic unloading in these elements. On the other hand, $$D_{eq}^p$$ increases monotonically for Elem. C, and reaches very high values when plastic instability occurs. Consequently, the plastic deformation localizes in this element. These trends may be correlated to the results obtained with the application of the Marciniak-Kuczynski approach (Ben Bettaieb and Abed-Meraim, 2015). Indeed, elements A and B play the role of the homogeneous zone, whereas element C plays the role of the imperfection band. In contrast to the Marciniak-Kuczynski analysis, strain localization can obviously be detected in the FE simulations without the need for introducing any initial geometric or material imperfection.

Figure 13
Evolution of the equivalent plastic strain rate as a function of the punch depth U: (a) Specimen with w = 30 mm; (b) Specimen with w = 90 mm; (c) Specimen with w = 185 mm.

One of the advantages of using a multiscale model at each integration point of the FE simulations (compared to phenomenological modeling) is the capability of accurately accounting for the effect of initial and induced plastic anisotropy during the deformation. This anisotropy has a significant influence on the mechanical behavior at the different integration points, and consequently on the overall response of the specimen. To further investigate this aspect, the final crystallographic textures at the different elements for the various specimens are presented in Figure 14, in the form of {111} pole figures. It is clearly revealed from this figure that the final texture of Elem. A (for all specimen geometries) is not significantly different from the initial texture shown in Figure 1 (b). This result is quite expectable considering the fact that the deformation in this element is the most moderate, as compared to the other analyzed elements (Elem. B and Elem. C). Therefore, the influence of deformation-induced anisotropy due to grain rotation during plastic deformation may not be quite important at the end of the simulation. On the other hand, the final textures corresponding to Elem. C are significantly different from the initial texture shown in Figure 1 (b). This result is also quite expectable considering the significant evolution of the strain rate in this element, as observed in Figure 13.

Figure 14
Final texture of the studied polycrystalline aggregate, {111} pole figure: (a) Elem. A, w = 30 mm ; (b) Elem. B, w = 30 mm; (c) Elem. C, w = 30 mm; (d) Elem. A, w = 90 mm ; (e) Elem. B, w = 90 mm; (f) Elem. C, w = 90 mm; (g) Elem. A, w = 185 mm ; (h) Elem. B, w = 185 mm; (i) Elem. C, w = 185 mm.

5 CONCLUDING REMARKS

In order to accurately simulate sheet forming processes, a suitable constitutive framework is required. To this aim, a multiscale polycrystalline model has been developed and implemented into the ABAQUS/Implicit finite element code. In this physically-based model, the link between the microscopic and macroscopic scales is achieved by using the self-consistent scale-transition scheme. The proposed model accounts for the initial and induced plastic anisotropy through the modeling of the texture evolution. The implemented constitutive model has been coupled with an indicator of loss of ellipticity in order to predict the occurrence of macroscopic plastic instability. This coupling has been used mainly in the following two applications:

• To investigate the effect of the number of grains, which constitute the polycrystalline aggregate, on the predicted forming limit diagram.

• To predict the onset and the evolution of plastic instability within sheet specimens undergoing a deep drawing process.

References

• Akpama, H. K., Ben Bettaieb, M., Abed-Meraim, F. (2016). Influence of the yield surface curvature on the forming limit diagrams predicted by crystal plasticity theory. Latin American Journal of Solids and Structures 13:2231-2250.
• Beaudoin, A. J., Dawson, P. R., Mathur, K. K., Kocks, U. F., Korzekwa, D. A. (1994). Application of polycrystal plasticity to sheet forming. Computer Methods in Applied Mechanics and Engineering 117:49-70.
• Ben Bettaieb, M., & Abed-Meraim, F. (2015). Investigation of localized necking in substrate-supported metal layers: comparison of bifurcation and imperfection analyses. International Journal of Plasticity 65:168-190.
• Ben Bettaieb, M., Débordes, O., Dogui, A., Duchêne, L., Keller, C. (2012). On the numerical integration of rate independent single crystal behavior at large strain. International Journal of Plasticity 32:184-217.
• Eshelby, J. D. 1957. The determination of the elastic field of an ellipsoidal inclusion, and related problems. In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. 241:376-396.
• Guan, Y., Pourboghrat, F., Barlat, F. (2006). Finite element modeling of tube hydroforming of polycrystalline aluminum alloy extrusions. International Journal of Plasticity 22:2366-2393.
• Heller, T., Engl, B., Ehrhardt, B., Esdohr, J. (1998). New high-strength steels-production, properties and applications. In Iron and Steel Society/AIME, 40th Mechanical Working and Steel Processing Conference Proceedings 36:25-34.
• Jung, K. H., Kim, D. K., Im, Y. T., Lee, Y. S. (2013). Prediction of the effects of hardening and texture heterogeneities by finite element analysis based on the Taylor model. International Journal of Plasticity 42:120-140.
• Kami, A., Dariani, B. M., Vanini, A. S., Comsa, D. S., Banabic, D. (2015). Numerical determination of the forming limit curves of anisotropic sheet metals using GTN damage model. Journal of Materials Processing Technology 216:472-483.
• Keeler, S. P., Backofen, W. A. (1963). Plastic instability and fracture in sheets stretched over rigid punches. Transactions of the ASM 56:25-48.
• Lipinski, P., Berveiller, M. (1989). Elastoplasticity of micro-inhomogeneous metals at large strains. International Journal of Plasticity 5:149-172.
• Lipinski, P., Berveiller, M., Reubrez, E., Morreale, J. (1995). Transition theories of elastic-plastic deformation of metallic polycrystals. Archive of Applied Mechanics 65:291-311.
• Mansouri, L. Z., Chalal, H., Abed-Meraim, F. (2014). Ductility limit prediction using a GTN damage model coupled with localization bifurcation analysis. Mechanics of Materials 76:64-92.
• Nakamachi, E. (1997). Optimum design of metal forming based on finite element analysis and nonlinear programming. In Tatsumi, T. et al. (Eds.), Theoretical and Applied Mechanics 1996:323-340.
• Nakamachi, E., Xie, C. L., Harimoto, M. (2001). Drawability assessment of BCC steel sheet by using elastic/crystalline viscoplastic finite element analyses. International Journal of Mechanical Sciences 43:631-652.
• Rudnicki, J. W., Rice, J. R. (1975). Conditions for the localization of deformation in pressure-sensitive dilatant materials. Journal of the Mechanics and Physics of Solids 23:371-394.
• Schmid, E., Boas, W. (1935). Kristallplastizität: mit besonderer Berücksichtigung der Metalle. Springer-Verlag 17.

Publication Dates

• Publication in this collection
  Oct 2017

History

• Received
  21 Nov 2016
• Accepted
  17 Mar 2017