Algorithmic Formulations of Evolutionary Anisotropic Plasticity Models Based on Non-Associated Flow Rule

Abstract

In the present paper, orthotropic elasto-plastic constitutive formulations for sheet metal forming based on non-associated flow rule that assume distortion of yield function/plastic potential with ongoing deformation process are analyzed. The yield function/plastic potential are considered as two different functions with functional form as orthotropic quadratic Hill or non-quadratic Karafillis-Boyce stress function. Based on the principle of plastic work equivalence, anisotropy parameters of the utilized yield function/plastic potential are set as functions of the equivalent plastic strain. In the constitutive formulation, for this internal variable, evolution equation consistent with the same principle of plastic work equivalence is introduced. For DC06 sheet sample with reported significant variation of the incremental r-values with straining, predictions of the evolution of the yield stress and r-value directional dependences with straining obtained by the analyzed models are presented. The algorithmic formulations of the analyzed constitutive models are derived by application of the implicit return mapping algorithm. For the derived stress integration procedures the accuracy is investigated by calculating iso-error maps. The maps are compared according to the flow rule and involved orthotropic stress functions. It has been revealed that although there is a difference in maps configuration there is no prominent difference in error magnitudes.

Keywords:
constitutive modeling; sheet metals; anisotropy evolution; non-associated flow rule; implicit return mapping

1 INTRODUCTION

Sheet metals are produced by complex thermo-mechanical processing routes among which the production rolling steps introduce a preferential orientation to the grains called texture. This preferential orientation results in anisotropic material behavior that is closely related to the formability of the material. The finite element codes that are widely used to design sheet forming processes most frequently utilize the phenomenological plasticity constitutive theories that are based on yield surface concept and plastic potential theory. The basis of these phenomenological constitutive models are: a yield criterion representing surface that separates the elastic and plastic regions in the stress space, a plastic potential whose gradient represents direction of the plastic strain rate and a hardening rule by which evolution of the yield surface is described.

The simplest and the most utilized approach in phenomenological plasticity theories that are intended for the metallic materials is the concept of isotropic hardening and use of the assumption known as associated flow rule or normality condition.

Associated flow rule assumption implies that the yield function and plastic potential function are identical and it is often regarded as the necessary condition for the metal plasticity fundamentals (stability of plastic flow and uniqueness). However, according to some theoretical observations (Mroz, 1963; Runesson and Mroz, 1989; Stoughton, 2002; Stoughton and Yoon, 2004, 2006; Cvitanić et al., 2008), an appropriate plastic flow description can be achieved by using less restrictive constrains over the constitutive equations than the associated flow rule. By non-associated flow rule, yielding and plastic flow are described by two different functions. In the last few years, several metal plasticity non-associated formulations were developed and this approach became particularly engaging in sheet metal plasticity. An acceptable constitutive model for a sheet material based on associated flow rule requires a yield function that can simultaneously predict anisotropy of yield stresses and anisotropy of plastic flow. Under the assumption of associated flow rule, various anisotropic phenomenological yield functions have been proposed for metallic sheets (Hill, 1948, 1979, 1990, 1993; Barlat and Lian, 1989; Karafillis and Boyce, 1993; Barlat et al., 1991, 1997, 2003, 2005; Banabic et al., 2003, 2005; Aretz, 2005; and others). Anisotropy is introduced into these stress functions by the parameters determined by the data related to the directional dependences of the yield stresses and Lankford parameters as well as data related to the biaxial tensile tests. Lankford parameter, also known as r-value or plastic strain ratio, is used as the measure of plastic flow. It is defined as the ratio of the sheet specimen transverse and thickness true plastic strain increments in uniaxial tensile testing. According to the standards, r-value is calculated by linear regression of the transverse versus longitudinal plastic strain plot between certain limits of the measured strains. Hence, it is considered as a constant value. In orthotropic sheet metal models based on non-associated flow rule, yield function and plastic potential function might have an identical functional form but their parameters can be defined by different experimental data. That means that the yield function can be adjusted to the yield stresses and the plastic potential can be adjusted to r-values. Such formulation can address sufficient number of experimental data and can result in acceptable predictions of the uniaxial material behavior even if yield/plastic potential stress functions with a low number of anisotropic parameters are utilized (Lademo et al., 1999; Stougthon, 2002; Stoughton and Yoon, 2004; Cvitanić et al., 2008; Safaei et al., 2014).

Isotropic hardening concept implies a proportional expansion of the yield surface with ongoing deformation process, without any change of its shape and position. In constitutive formulations intended for sheet metals, by application of isotropic hardening and standard calculation procedure for r-value, anisotropy parameters are calculated using the initial yield stress ratios and/or constant r-values. Such approach results in fixed anisotropy parameters and accordingly any possible distortion of yield function/plastic potential with continuation of the plastic deformation process is neglected. In addition, uniform size change of the orthotropic yield function in the stress space is controlled by single scalar hardening variable and mostly flow stress curve for the rolling direction is adopted as the expansion rule. However, some recent experimental studies report alternation of the yield stress ratios and/or instantaneous r-values with evolution of sheet texture during deformation process (Zamiri and Pourboghrat, 2007; An et al., 2013; Safaei et al., 2014). Furthermore, numerous studies related to the application of the orthotropic plasticity formulations with constant anisotropy parameters in predicting complex forming processes indicate that possible model improvements could be achieved by incorporating the evolution of yield stress ratios and r-values into the model. This is particularly evident in the simulations of the cylindrical cup drawing problem. In this forming process, the plastic anisotropy has important effect on the formation of undulative top edge of the finally drawn cup often termed as earing. Numerous phenomenological orthotropic yield functions with constant (fixed) anisotropy parameters under associated or non-associated flow rule have been tested in this forming process (Yoon et al., 1998, 1999, 2000, 2004, 2006; Cvitanić et al., 2008; Taherizadeh et al., 2010, 2011; Park and Chung, 2012; Safaei et al., 2013; Vrh et al., 2014; and others). These extensive numerical studies clearly indicate that the predicted earing trend (location of peaks and valleys) is the mirror image of the predicted r-value trend with the respect to the transverse direction, and that the prediction of the earing profile amplitude (the maximum difference in cup heights) is correlated to the prediction of the yield stress anisotropy amplitude. Considering predictions of the complex phenomenological yield functions that are capable to accurately predict directional dependence of the r-values, it can be observed that there is still discrepancy between their predictions of earing profile and experimental cup heights. Although the description of the material behavior is not the only factor that influences the finite element predictions, due to the expressed correlation between predicted r-values and predicted earing profile, it can be concluded that reliable material model should take into account not only directional dependences of the uniaxial material data determined at the initial yield state or their averaged values but also their possible evolution with continuation of the deformation process. That means that the anisotropy parameters of the yield function/plastic potential should be altered by a certain measure of plastic flow.

In this paper, constitutive formulations based on non-associated flow rule and distortional evolution of the orthotropic yield function/plastic potential with ongoing deformation process are considered. As yield function/plastic potential orthotropic four parametric quadratic stress function proposed by Hill (1948) or non-quadratic stress function proposed by Karafillis-Boyce (1993) are utilized. Anisotropy parameters of the utilized yield function/plastic potential are set as polynomial functions of the equivalent plastic strain and evolution equation for this hardening parameter consistent with the principle of the plastic work equivalence is introduced. Previously developed algorithmic formulations of the evolutionary anisotropic model based on the considered stress functions and associated flow rule (Cvitanić et al., 2016) are now extended for the case of non-associated flow rule. Considered algorithmic formulations are derived by application of the implicit return mapping algorithm. Numerical analysis of accuracy of the proposed stress integration procedures is performed and compared according to the flow rule and the functional form of the involved stress functions by calculating iso-error maps. The constitutive descriptions are developed and analyzed considering experimental data for DC06 steel sheet reported by Safaei et al. (2014).

2 EVOLUTIONARY ANISOTROPIC PLASTICITY FORMULATION BASED ON NON-ASSOCIATED FLOW RULE

Assuming isotropic linear elasticity and additive decomposition of the strain tensor increment d( into elastic _^d(e and plastic part _^d(p , the stress tensor increment dσ reads

$$d\sigma =C^e:d{\epsilon }^e=C^e:(d\epsilon -d{\epsilon }^p)$$ (1)

where C^e is the tensor of elastic module. Considering sheet material with anisotropy evolution, the yield criterion is stated as follows

$$F=f_y(\sigma ,{\overline{\epsilon }}^p)-\kappa ({\overline{\epsilon }}^p)=0$$ (2)

where $$f_y(\sigma ,{\overline{\epsilon }}^p)$$ is an orthotropic yield function with orthotropic parameters introduced as functions of hardening parameter $${\overline{\epsilon }}^p$$ and $$\kappa ({\overline{\epsilon }}^p)$$ is a scalar function representing stress-strain relation for the referent direction. According to the plastic potential theory, the plastic part of the strain tensor increment d ε ^p is proportional to the gradient of the stress function named plastic potential function

$$d{\epsilon }^p=d\lambda \frac{\partial f_p(\sigma ,{\overline{\epsilon }}^p)}{\partial \sigma }$$ (3)

where dλ is a non-negative scalar called plastic multiplier or consistency parameter. Considering anisotropy evolution, plastic potential $$f_p(\sigma ,{\overline{\epsilon }}^p)$$ is also introduced as an orthotropic stress function with orthotropic parameters stated as functions of the parameter $${\overline{\epsilon }}^p.$$ If the plastic potential and yield function are identical $$f_p(\sigma ,{\overline{\epsilon }}^p)\equiv f_y(\sigma ,{\overline{\epsilon }}^p)$$, yielding and plastic flow are described by the same function and Eq. (3) becomes the so-called associated flow rule. In the present formulation, the parameter $${\overline{\epsilon }}^p$$ is considered as equivalent plastic strain that obeys the following form of the principle of plastic work equivalence

$$f_y(\sigma ,{\overline{\epsilon }}^p)d{\overline{\epsilon }}^p=\sigma :d{\epsilon }^p$$ (4)

If the plastic potential function fulfils identity $$f_p(\sigma ,{\overline{\epsilon }}^p)=\sigma :\partial f_p(\sigma ,{\overline{\epsilon }}^p)/\partial \sigma$$, using Eqs. (3) and (4), following evolution equation for the parameter $${\overline{\epsilon }}^p$$ is obtained

$d{\overline{\epsilon }}^p=d\lambda \frac{\sigma :\partial f_p(\sigma ,{\overline{\epsilon }}^p)/\partial \sigma }{f_y(\sigma ,{\overline{\epsilon }}^p)}=d\lambda \frac{f_p(\sigma ,{\overline{\epsilon }}^p)}{f_y(\sigma ,{\overline{\epsilon }}^p)}$ (5)

where for the associated flow rule $$d{\overline{\epsilon }}^p=d\lambda$$. If deformation process is elastic, the incremental changes of the internal variables (plastic strain tensor d( ^P and parameter $${\overline{\epsilon }}^p$$) vanish and dλ =0. Therefore, for the hardening material the plastic multiplier obeys the complementary conditions dλ ≥ 0, F ≤ 0, dλF = 0 and consistency condition dλdF = 0.

3 ANISOTROPIC STRESS FUNCTIONS WITH ANISOTROPY EVOLUTION

The analyzed constitutive formulations utilize yield function and plastic potential with functional form of orthotropic quadratic Hill (1948) stress function or non-quadratic Karafillis-Boyce (1993) stress function. In sheet metal forming it is common practice to assume that the sheet is approximately subjected to plane stress conditions and that material exhibits orthotropic symmetry in plastic properties. Accordingly, the utilized yield functions/plastic potentials are stated as functions of in-plane stress components σ_xx , σ_yy and σ_zz , where x-axis denotes the original sheet rolling direction and y-axis denotes the direction in sheet plane transverse to the rolling direction. The z-axis denotes the sheet normal direction.

The utilized Hill stress function and Karafillis-Boyce stress function inclose four parameters by which material plastic anisotropy can be described. Under the associated flow rule, the utilized stress function acts as yield function as well as plastic potential, therefore, its anisotropy parameters can be adjusted to the yield stresses or experimental data indicating plastic flow. In the non-associated formulation, anisotropy parameters of the Hill or Karafillis-Boyce yield function are defined in terms of three directional yield stresses obtained in the uniaxial tension of the specimens oriented at 0(, 45( and 90( to the rolling direction and equibiaxial yield stress. The associated yield stresses are denoted as σ₀, σ₄₅, σ₉₀ and σ_b . The anisotropy parameters of the Hill or Karafillis-Boyce plastic potential are defined in terms of experimental data indicating plastic flow such as Lankford parameter that reads

$$r=\frac{d{\epsilon }_{22}^p}{d{\epsilon }_{33}^p}$$ (6)

where $$d{\epsilon }_{22}^p$$ and $$d{\epsilon }_{33}^p$$ are width and thickness plastic strain increments, respectively, obtained in uniaxial sheet specimen tension. In calculating parameters of the plastic potential, three plastic strain ratios obtained in the uniaxial tensions along 0(, 45( and 90( to the rolling direction and the yield stress for the referent direction σ₀ are used. The associated plastic strain ratios are denoted as r ₀, r ₄₅ and r ₉₀.

3.1 Hill Stress Function and Karafillis-Boyce Stress Function

The orthotropic quadratic Hill (1948) stress function is derived as an extension of the isotropic von Mises yield function. For plane stress condition Hill stress function can be written in the following form

$f_y=\sqrt{{\lambda }_1{\sigma }_{xx}^2+{\lambda }_2{\sigma }_{yy}^2-2\nu {\sigma }_{xx}^{}{\sigma }_{yy}^{}+2\rho {\sigma }_{xy}^2}={\sigma }_y^{}$ (7)

where λ₁, λ₂, v and ρ are anisotropic material parameters that can be adjusted to experimental data and σ_y is the yield stress for the referent direction. For Hill stress function calculation procedure results in the explicit expressions for anisotropy parameters. The parameters of the Hill stress function adjusted to the yield stresses read

$$\begin{array}{l}{\lambda }_1=\mathrm{1,}{\lambda }_2=({\sigma }_0/{\sigma }_{90}{)}^2,\nu =0.5(1+({\sigma }_0/{\sigma }_{90}{)}^2+({\sigma }_0/{\sigma }_b{)}^2),\\ \rho =2({\sigma }_0/{\sigma }_{45}{)}^2-0.5({\sigma }_0/{\sigma }_b{)}^2\end{array}$$ (8)

If the Hill stress function is adjusted to the plastic strain ratios, the calculation procedure results in following explicit expressions

$${\lambda }_1=\mathrm{1,}{\lambda }_2=\frac{1+1/r_{90}}{1+1/r_0},\nu =\frac{1}{1+1/r_0},\rho =\frac{(1+2r_{45})(1/r_0+1/r_{90})}{2(1+1/r_0)}$$ (9)

The orthotropic Karafillis-Boyce (1993) stress function is a linear combination of two convex non-quadratic functions

$$f_y={\left(\frac{\left(1-c\right)}{2}\left({|{\tilde{s}}_1-{\tilde{s}}_2|}^m+{|{\tilde{s}}_2-{\tilde{s}}_3|}^m+{|{\tilde{s}}_3-{\tilde{s}}_1|}^m\right)+\frac{c}{2^{}}\frac{3^m}{\left(2^{m-1}+1\right)}\left({|{\tilde{s}}_1|}^m+|{\tilde{s}}_2|{}^m+{|{\tilde{s}}_3|}^m\right)\right)}^{1/m}={\sigma }_y^{}$$ (10)

where c is a weighting parameter and exponent m is an even number. In the above equation, ${\tilde{s}}_1$, ${\tilde{s}}_2$ and ${\tilde{s}}_3$ are the principal values of the so-called isotropic plasticity equivalent stress tensor. For plane stress conditions, these values can be calculated as

${\tilde{s}}_{\mathrm{1,2}}=\frac{{\tilde{s}}_{xx}+{\tilde{s}}_{yy}}{2}\pm \sqrt{{\left(\frac{{\tilde{s}}_{xx}-{\tilde{s}}_{yy}}{2}\right)}^2+{\tilde{s}}_{xy}^2},{\tilde{s}}_3={\tilde{s}}_{zz}$ (11)

where

$\left\{\begin{array}{c}{\tilde{s}}_{xx}\\ {\tilde{s}}_{yy}\\ {\tilde{s}}_{zz}\\ {\tilde{s}}_{xy}\end{array}\right\}=C\left[\begin{array}{cc}1& \begin{array}{ccc}{\beta }_1& {\beta }_2& 0\end{array}\\ \begin{array}{c}{\beta }_1\\ \begin{array}{c}{\beta }_2\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{ccc}{\alpha }_1& {\beta }_3& 0\end{array}\\ \begin{array}{c}\begin{array}{ccc}{\beta }_3& {\alpha }_2& 0\end{array}\\ \begin{array}{ccc}0& 0& {\gamma }_3\end{array}\end{array}\end{array}\end{array}\right]\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{zz}=0\\ {\sigma }_{xy}\end{array}\right\}$ (12)

and

${\beta }_1=\frac{{\alpha }_2-{\alpha }_1-1}{2},{\beta }_2=\frac{{\alpha }_1-{\alpha }_2-1}{2},{\beta }_3=\frac{1-{\alpha }_1-{\alpha }_2}{2}$ (13)

In Karafillis-Boyce stress function material anisotropy is defined by the parameters C, α ₁, α ₂ and γ ₃. For isotropic material, these parameters have the values C = 2/3, α ₁ = α ₂ =1 and γ ₃ = 3/2 and stress components defined by Eq. (12) reduce to the components of the stress deviator tensor.

For Karafillis-Boyce stress function, calibration procedure i.e. calculation of anisotropy parameters, leads to the system of non-linear equations that can be solved using a numerical iterative procedure.

3.2 Evolution of Hill and Karafillis-Boyce Yield Function/Plastic Potential

In this paper, experimentally determined directional dependences of the uniaxial material properties and their evolution with ongoing plastic deformation process for 0.8 mm thick DC06 steel sheet reported by Safaei et al. (2014) are used to develop and analyze evolutionary anisotropic constitutive models. The utilized data are related to the uniaxial straining tests of the seven specimens with orientations 0°, 15°, 30°, 45°, 60°, 75° and 90° to the rolling direction. The reported data utilized in model development are: 1) parameters of the combined Swift-Voce hardening law by which the experimental true stress and longitudinal true plastic strain are approximated; 2) parameters of the 3^rd order polynomial fit by which the experimental transverse and longitudinal true plastic strains are approximated. In the previous study conducted by the authors (Cvitanić et al., 2016), using these data, parameters of the Hill and Karafillis-Boyce yield function/plastic potential are introduced as polynomial functions of the equivalent plastic strain. In the first step, following procedure based on the principle of plastic work equivalence presented in Safaei et al. (2014), for each specimen orientation yield stresses σ_y and instantaneous r-values corresponding to the certain amount of the equivalent true plastic strain are calculated. In the adopted approach, longitudinal true plastic strain in the rolling direction is used as the equivalent true plastic strain. In calculating instantaneous r-values incompressibility hypothesis is applied, thus Eq. (6) is considered in the following form

$$r_{\theta }=-\frac{d{\epsilon }_{\theta +90}^p}{d{\epsilon }_{\theta +90}^p+d{\epsilon }_{\theta }^p}=-\frac{m_{\theta }}{m_{\theta }+1}$$ (14)

where $d{\epsilon }_{\theta }^p$ and $$d{\epsilon }_{\theta +90}^p$$ are the increments of true longitudinal and transverse plastic strains corresponding to the loading direction θ and direction θ + 90°, respectively. According to Eq. (14), instantaneous r-values corresponding to the certain longitudinal true plastic strain can be calculated using the slope _^mθ of the available 3^rd order polynomial fit by which the experimental transverse and longitudinal true plastic strains are approximated.

According to the adopted calculation procedure, yield stresses and r-values corresponding to the seven orientations (0°, 15°, 30°, 45°, 60°, 75° and 90°) and amounts of the equivalent plastic strain starting from 0.001 to 0.301 at each 0.002 increment are calculated. The calculated yield stress ratios (yield stresses normalized with yield stress for the rolling direction) and r-values corresponding to the selected values of the equivalent plastic strain are presented in Figure 1 and Figure 2, respectively. From Figure 1 it can be observed that the directional dependence trend of the yield stress ratios at the start of plastic deformation process ($${\overline{\epsilon }}^p=0.001$$) is rather distorted with further straining. Furthermore, as presented in Figure 2, for the considered material, there is a significant decrease of the r-values and certain alternation of the r-value directional dependence with ongoing deformation.

Figure 1
Yield stress ratio directional dependences corresponding to several values of the equivalent plastic strain. Predictions of yield stress ratios obtained by Hill yield function.

Figure 2
Lankford parameter directional dependences corresponding to several values of the equivalent plastic strain. Predictions of Lankford parameters obtained by Hill plastic potential.

These observations clearly indicate the suitability of the non-associated flow rule approach, by which the evolution of yield stress and r-value directional dependences can be independently considered and described. In Figures 1 and 2, predictions of the yield stress and r-value directional dependences corresponding to the selected equivalent plastic strains obtained by the Hill yield function/plastic potential are also presented. The predictions of the Karafillis-Boyce yield function/plastic potential are almost identical to the predictions obtained by the Hill functions and therefore they are not separately presented.

The anisotropy parameters of Hill and Karafillis-Boyce yield functions are calculated using the yield stresses corresponding to the orientations 0(, 45( and 90( (σ_0, σ₄₅, σ₉₀) and assuming that the yield stress at balanced biaxial stress state is the averaged value of the yield stresses corresponding to the rolling and transverse direction σ_{b =} (σ_{0 +} σ₉₀ ) / 2. The parameters of the Hill and Karafillis-Boyce plastic potentials are calculated using r-values corresponding to the orientations 0(, 45( and 90( (r ₀, r ₄₅, r ₉₀). The parameters of the analyzed Hill and Karafillis-Boyce yield function/plastic potential are calculated considering calculated yield stresses and r-values corresponding to the equivalent plastic strain values starting from 0.001 to 0.301 at each 0.002 increment. From Figures 1 and 2, it can be observed that the analyzed yield functions poorly predict the pronounced directional dependence at the onset of the deformation process while the analyzed plastic potentials result in good predictions of the r-value directional dependences, particularly for greater strain levels.

In Figure 3 predictions of yield stress ratios obtained by the models based on the associated flow rule and Hill or Karafillis-Boyce stress function are presented. In these models, utilized orthotropic stress function adjusted to r-values acts as plastic potential and as yield function. It can be observed that the analyzed associated models fail in predicting directional dependence of the yield stress ratios for the considered deformation levels. These results in addition support the use of the non-associated flow rule approach coupled with analyzed four parametric stress functions in describing plastic anisotropy and its evolution for the considered material.

In order to relate the anisotropy parameters with the equivalent plastic strain, the fourth order polynomial fit is used

$$Poly4\left({\overline{\epsilon }}^p\right)=a+b_1\cdot {\overline{\epsilon }}^p+b_2{\left({\overline{\epsilon }}^p\right)}^2+b_3{\left({\overline{\epsilon }}^p\right)}^3+b_4{\left({\overline{\epsilon }}^p\right)}^4$$ (15)

Figure 3
Yield stress ratio directional dependences predicted by associated model and a) Hill ; b) Karafillis-Boyce function adjusted to r-values corresponding to several values of the equivalent plastic strain.

The calculated polynomial parameters for the analyzed yield functions/plastic potentials are obtained using the least square method and are presented in Table 1. Figures 4 and 5 show values of each anisotropy parameter corresponding to the selected values of the equivalent plastic strain and related polynomial fit. From these figures it can be seen that there is a good correlation between adopted fit and plastic potential parameters, while there is a certain discrepancy for yield function parameters for the lower plastic levels.

Figure 4
Poly4 function fit for the anisotropy parameters of a) Hill; b) Karafillis-Boyce yield function.

Figure 5
Poly4 function fit for the anisotropy parameters of a) Hill; b) Karafillis-Boyce plastic potential.

Table 1
Parameters of Poly4 fit for Hill and Karafillis-Boyce yield function/plastic potential.

	Hill yield function				Hill plastic potential function
	λ1	λ2	ν	ρ	λ1	λ2	ν	ρ
a	1	0.92024	0.481464	1.429272	1	0.929330	0.67630	1.518681
b 1	0	3.52471	0.827483	4.650582	0	-0.34531	-0.26861	0.211961
b 2	0	-47.34067	-11.21203	-66.5226	0	0.455731	-0.05202	-0.5717
b 3	0	225.4881	53.614063	323.6818	0	-0.15724	-0.00091	-1.00629
b 4	0	-351.6768	-83.80070	-512.050	0	0.219299	0.00151	1.586853
	Karafillis-Boyce yield function				Karafillis-Boyce plastic potential function
	C	α1	α2	γ3	C	α1	α2	γ3
a	0.6664	0.95801	0.9768	1.4632	0.6631	0.97304	0.854890	1.511186
b 1	0.0074	1.84314	1.0031	2.3715	0.0108	-0.12946	0.055896	0.016846
b 2	-0.0871	-24.6067	-13.2496	-33.751	-0.004	0.19650	0.245121	-0.13490
b 3	0.3861	116.8851	62.629	163.903	-0.009	-0.03297	-0.08198	-0.17842
b 4	-0.5770	-182.018	-97.2598	-259.026	0.0066	0.01965	0.038359	0.315813

Contours of the analyzed yield functions and plastic potentials in normalized stress space for zero shear stress corresponding to the selected values of the equivalent plastic strain are presented in Figure 6 and Figure 7, respectively. For material that obeys the isotropic hardening concept, yield contours that correspond to the different amount of plastic strain should coincide if presented in normalized stress space. Considering obtained contours, it can be observed that there is a significant discrepancy between initial and subsequent yield contours and notable evolution of plastic potential contour particularly for the quadratic potential function.

IMG35 $${\overline{\epsilon }}^p=0.001;0.101;0.201;0.301$$

Figure 6
Contours of a) Hill; b) Karafillis-Boyce yield function corresponding to IMG35.

IMG36 $${\overline{\epsilon }}^p=0.001;0.101;0.201;0.301$$

Figure 7
Contours of a) Hill; b) Karafillis-Boyce plastic potential corresponding to IMG36.

4 STRESS INTEGRATION PROCEDURE

In this section, for the analyzed evolutionary anisotropic plasticity models based on non-associated flow rule, computational procedure for calculating state variables at time _{^tn +1} ($${\sigma }_{n+1},{\overline{\epsilon }}_{n+1}^p$$) based on the known state variables at time _^tn ($${\sigma }_n,{\overline{\epsilon }}_n^p$$) and known increment of total deformation Δε is derived. The derived procedure is based on the implicit return mapping (Simo and Hughes, 1988; Yoon et al., 1999; Cvitanić et al., 2008) and presents an extension of the procedure previously derived for evolutionary anisotropic model based on the associated flow rule (Cvitanić et al., 2016).

By the application of the implicit return mapping procedure, the stress solution is obtained in two steps. In the elastic predictor step, the strain increment ∆( is assumed to be elastic and trial elastic stress tensor _^σtrial is calculated based on the previously converged solution σ_n

$${\sigma }^{trial}={\sigma }_n+C^e:\Delta \epsilon$$ (16)

If the trial state violates the yield condition, the plastic correction step is performed assuming trial state as initial condition. In this step the final stress _{^σn +1} is stated as

$${\sigma }_{n+1}={\sigma }_n+C^e:(\Delta \epsilon -\Delta {\epsilon }^p)={\sigma }^{trial}-C^e:\Delta {\epsilon }^p$$ (17)

and evolution equations for the plastic strain tensor and hardening parameter are integrated to restore the consistency condition. By application of implicit Euler backward integration procedure and assuming non-associated flow rule, increment of plastic strain tensor is approximated as follows

$\Delta {\epsilon }^p=\Delta \lambda {(\partial f_p(\sigma ,{\overline{\epsilon }}^p)/\partial \sigma )}_{n+1}$ (18)

Consistently, the increment of the hardening parameter defined by Eq. (5) is approximated as follows

$\Delta {\overline{\epsilon }}^p=\Delta \lambda {(f_p(\sigma ,{\overline{\epsilon }}^p)/f_y(\sigma ,{\overline{\epsilon }}^p))}_{n+1}$ (19)

In above approximations, ∆λ is an incremental consistency parameter that obeys discrete form of the complementary conditions

$\Delta \lambda \ge \mathrm{0,}F({\sigma }_{n+1},{\overline{\epsilon }}_{n+1}^p)\le \mathrm{0,}\Delta \lambda F({\sigma }_{n+1},{\overline{\epsilon }}_{n+1}^p)=0$ (20)

By using Eqs. (16)-(19), incremental form of the constitutive model can be stated by the following system of four non-linear algebraic equations (tensorial and scalar)

$${\Phi }_1\equiv f_y({\sigma }_{n+1},{\overline{\epsilon }}_n^p+\Delta {\overline{\epsilon }}^p)-{\rho }_{n+1}=0$$ (21)

${\Phi }_2\equiv {(C^e)}^{-1}:({\sigma }_{n+1}-{\sigma }^{trial})+\Delta \lambda n_{n+1}=0$ (22)

${\Phi }_3={\rho }_{n+1}-\kappa ({\overline{\epsilon }}_n^p+\Delta {\overline{\epsilon }}^p)=0$ (23)

$${\Phi }_4=f_y({\sigma }_{n+1},{\overline{\epsilon }}_n^p+\Delta {\overline{\epsilon }}^p)\Delta {\overline{\epsilon }}^p-f_p({\sigma }_{n+1},{\overline{\epsilon }}_n^p+\Delta {\overline{\epsilon }}^p)\Delta \lambda =0$$ (24)

where $n=\partial f_p(\sigma ,{\overline{\epsilon }}^p)/\partial \sigma$. The solution of above equation system (four unknowns $${\sigma }_{n+1},{\rho }_{n+1},\Delta {\overline{\epsilon }}^p,\Delta \lambda$$) can be obtained using iterative Newton-Raphson procedure. At each iteration, denoted as k, above four equations are linearized around the current values of state variables to obtain

$${\Phi }_1^{(k)}+m^{(k)}:\Delta {\sigma }^{(k)}+{(\partial f_y/\partial {\overline{\epsilon }}^p)}^{(k)}\delta \Delta {\overline{\epsilon }}^{p(k)}-\Delta {\rho }^{(k)}=0$$ (25)

$${\Phi }_2^{(k)}+{(C^e)}^{-1}:\Delta {\sigma }^{(k)}+\Delta {\lambda }^{(k)}{(\partial n/\partial \sigma )}^{(k)}:\Delta {\sigma }^{(k)}+\Delta {\lambda }^{(k)}{(\partial n/\partial {\overline{\epsilon }}^p)}^{(k)}\delta \Delta {\overline{\epsilon }}^{p(k)}+\delta \Delta {\lambda }^{(k)}{n}^{(k)}=0$$ (26)

$${\Phi }_3{}^{(k)}+\Delta {\rho }^{(k)}-{(d\kappa /d{\overline{\epsilon }}^p)}^{(k)}\delta \Delta {\overline{\epsilon }}^{p(k)}=0$$ (27)

$\begin{array}{l}{\Phi }_4^{\left(k\right)}+\Delta {\overline{\epsilon }}^{p\left(k\right)}{m}^{\left(k\right)}:\Delta {\sigma }^{\left(k\right)}+\Delta {\overline{\epsilon }}^{p\left(k\right)}{\left(\partial f_y/\partial {\overline{\epsilon }}^p\right)}^{\left(k\right)}\delta \Delta {\overline{\epsilon }}^{p\left(k\right)}+j_y^{\left(k\right)}\delta \Delta {\overline{\epsilon }}^{p\left(k\right)}\text{}\\ -\Delta {\lambda }^{\left(k\right)}{n}^{\left(k\right)}:\Delta {\sigma }^{\left(k\right)}-\Delta {\lambda }^{\left(k\right)}{\left(\partial f_p/\partial {\overline{\epsilon }}^p\right)}^{\left(k\right)}\delta \Delta {\overline{\epsilon }}^{p\left(k\right)}-f_p^{\left(k\right)}\delta \Delta {\lambda }^{\left(k\right)}=0\end{array}$ (28)

where $m=\partial f_y(\sigma ,{\overline{\epsilon }}^p)/\partial \sigma$. By solving above linearized system, explicit expressions for the increments of state variables are obtained

$$\delta \Delta {\epsilon }^{p(k)}=\frac{{\Phi }_1^{(k)}-m^{(k)}:{(M^{(k)})}^{-1}:{\Phi }_2^{(k)}+{\Phi }_3^{(k)}-m^{(k)}:{(M^{(k)})}^{-1}:n^{(k)}({\Phi }_4^{(k)}/f_p^{(k)})}{\left({\left(\frac{d\kappa }{d{\overline{\epsilon }}^p}\right)}^{(k)}-{\left(\frac{\partial f_y}{\partial {\overline{\epsilon }}^p}\right)}^{(k)}\right)+m^{(k)}:{(M^{(k)})}^{-1}:n^{(k)}{R}^{(k)}+m^{(k)}:{(M^{(k)})}^{-1}:{\left(\frac{\partial n}{\partial {\overline{\epsilon }}^p}\right)}^{(k)}\Delta {\lambda }^{(k)}}$$ (29)

$$\Delta {\sigma }^{(k)}=-{(M^{(k)})}^{-1}:(n^{(k)}({\Phi }_4^{(k)}/f_p^{(k)})+{\Phi }_2^{(k)}+(n^{(k)}{R}^{(k)}+\Delta {\lambda }^{(k)}{(\partial n/\partial {\overline{\epsilon }}^p)}^{(k)})\delta \Delta {\overline{\epsilon }}^{p(k)})$$ (30)

$\delta \Delta {\lambda }^{(k)}=(1/f_p^{(k)})G^{(k)}:\Delta {\sigma }^{(k)}+({\Phi }_4^{(k)}/f_p^{(k)})+R^{(k)}\delta \Delta {\overline{\epsilon }}^{p(k)}$ (31)

$\Delta {\rho }^{(k)}=-{\Phi }_3^{(k)}+{(d\kappa /d{\overline{\epsilon }}^p)}^{(k)}\delta \Delta {\overline{\epsilon }}^{p(k)}$ (32)

where

$\begin{array}{l}{\overline{C}}^{\left(k\right)}={\left(C^e\right)}^{-1}+\Delta {\lambda }^{\left(k\right)}{\left(\partial n/\partial \sigma \right)}^{\left(k\right)},G^{\left(k\right)}=\Delta {\overline{\epsilon }}^{p\left(k\right)}{m}^{\left(k\right)}-\Delta {\lambda }^{\left(k\right)}{n}^{\left(k\right)},\text{}\\ M^{\left(k\right)}=n^{\left(k\right)}\otimes G^{\left(k\right)}/f_p^{\left(k\right)}+{\overline{C}}^{\left(k\right)},\text{}\\ R^{\left(k\right)}=\left(f_y^{\left(k\right)}/f_p^{\left(k\right)}\right)+\left(\Delta {\overline{\epsilon }}^{p\left(k\right)}{\left(\partial f_y/\partial {\overline{\epsilon }}^p\right)}^{\left(k\right)}-\Delta {\lambda }^{\left(k\right)}{\left(\partial f_p/\partial {\overline{\epsilon }}^p\right)}^{\left(k\right)}\right)/f_p^{\left(k\right)}.\text{}\end{array}$

Finally, updated state variables for non-associated formulation are defined as

$${\left\{\begin{array}{c}\Delta {\overline{\epsilon }}^p\\ \rho \\ \sigma \\ \Delta \lambda \end{array}\right\}}_{(n+1)}^{(k+1)}={\left\{\begin{array}{c}\Delta {\overline{\epsilon }}^p\\ \rho \\ \sigma \\ \Delta \lambda \end{array}\right\}}_{(n+1)}^{(k)}+{\left\{\begin{array}{c}\delta \Delta {\overline{\epsilon }}^p\\ \Delta \rho \\ \Delta \sigma \\ \delta \Delta \lambda \end{array}\right\}}_{(n+1)}^{(k)}$$ (33)

If the associated flow rule is assumed, according to Eq. (19) $\Delta {\overline{\epsilon }}^p=\Delta \lambda$ holds and integration scheme proposed by Cvitanić et al. (2016) is reproduced. On the other hand, if anisotropy evolution is neglected orthotropic parameters of yield function/plastic potential are not altered by hardening parameter and integration scheme conform the one developed by Cvitanić et al. (2008).

5 NUMERICAL ANALYSIS OF ACCURACY. ISO-ERROR MAPS

In order to estimate the accuracy of the proposed algorithm based on non-associated flow rule and distortional evolution of the yield function and the plastic potential function, iso-error maps (Simo and Hughes, 1988) are calculated. The iso-error maps are generally accepted as an effective and reliable tool for assessing the accuracy of the constitutive integration algorithms under realistic strain steps. For two-dimensional implementation of the developed algorithms, iso-error maps are generated following known procedure. For the selected stress point on the yield surface sequence of two- component strain increments is applied and for each strain increment numerical solution is obtained according to the tested algorithm. The solution obtained by the same algorithm by dividing the considered strain increment into large number of sub-increments of the same size is adopted as exact solution. Furthermore, for each considered strain increment, based on the obtained numerical and exact solution, the error is calculated as the percentage of the relative root mean square of errors between the computed stress tensor σ and exact stress tensor σ*. This error measure is defined as

$$\delta (\%)=\frac{\sqrt{(\sigma -{\sigma }^{\ast }):(\sigma -{\sigma }^{\ast })}}{\sqrt{{\sigma }^{\ast }:{\sigma }^{\ast }}}\cdot 100$$ (34)

Finally, for the considered stress point iso-error map is drawn as the contour plot of the error field presented over the two-component strain space. By inspection of these plots obtained for the representative stress points, dependence of the algorithm accuracy on the strain increment size can be assessed.

For the analyzed formulations, iso-error maps are calculated at three representative stress points on the yield surface: A-uniaxial, B-balanced biaxial and C- pure shear as shown in Figure 8.

Figure 8
Plane stress yield surface and points A, B and C for iso-error maps.

The strain increments along sheet rolling and transverse direction ranging from zero to six times of the yield strain _^(y are applied to the considered stress points. The calculations are performed assuming _^∆(xy = 0 thus _^σxy = 0 applies. As the exact solution, the stresses obtained by 100 sub-steps of each strain increment are utilized for the Hill formulations. For the Karafillis-Boyce formulations solutions obtained by 50 sub-steps are adopted.

In the tested algorithms, the polynomial relations for the anisotropic parameters of the Hill stress functions and Karafillis-Boyce stress functions corresponding to DC06 steel sheet that are presented in Section 3 are utilized. Furthermore, combined Swift-Voce stress-strain relation for the rolling direction reported by Safaei et al. (2014) for the considered material is used to govern the yield surface expansion

$$\kappa ({\overline{\epsilon }}^p)=c^{\prime }(k{({\overline{\epsilon }}_0^p+{\overline{\epsilon }}^p)}^n)+(1-c^{\prime })(R+Q(1-e^{-b{\overline{\epsilon }}^p}))$$ (35)

where c' = 0.848; k = 539.542 img = 0.012; n = 0.326; R = 29.247; Q = 557.223; b = 34.822. For the elastic constants, Young’s modulus E = 200 GPa and Poisson’s coefficient ν = 0.3 are utilized.

Figures 9 and 10 show calculated iso-error maps for three representative stress points obtained by the Hill and Karafillis-Boyce formulations, respectively. Calculated iso-error maps are compared to the maps obtained by the algorithm based on associated flow rule and distortional hardening concept coupled with the yield function adjusted to r-values (Cvitanić et al., 2016). For the analyzed formulations, the magnitudes of errors are comparable to the results reported in literature obtained by the associated/non-associated formulations based on the utilized orthotropic stress functions and isotropic hardening concept (Cvitanić et al. 2008). Commonly, considering points B and C, for the isotropic formulations the exact solutions are obtained for the loading directed along the yield surface symmetry axis. As shown in Figures 9 and 10, considering these points, for the analyzed orthotropic formulations the axes of exact solution are shifted to the stress symmetry axes. From the comparison of the iso-error maps for different stress points, it can be seen that the errors are relatively smaller for the balanced biaxial stress state for all analyzed formulations. Considering different flow rule approach, for the Hill formulations there is no significant difference between iso-error maps obtained by the associated or non-associated flow rule. For the Karafillis-Boyce formulations, the errors are relatively smaller if associated flow rule is utilized. Furthermore, from the comparison of the iso-error maps for different orthotropic stress functions it can be seen that although there is a difference in maps configuration there is no prominent difference in error magnitudes. Finally, it can be stated that for the analyzed formulations, a reduction in error magnitude can be expected when the strain size increment is reduced.

Figure 9
Iso-error maps based on Hill associated formulation for points: a) A; b) B; c) C and on Hill non-associated formulation for points: d) A; e) B; f) C.

Figure 10
Iso-error maps based on Karafillis-Boyce associated formulation for points: a) A; b) B; c) C and on Karafillis-Boyce non-associated formulation for points: d) A; e) B; f) C.

6 CONCLUSIONS

In the present paper, constitutive formulations for sheet metal forming based on non-associated flow rule that enable distortion of the yield function/plastic potential with on-going deformation process are analyzed. The formulations are developed considering evolution of the yield stress ratios and r-values with straining obtained in uniaxial tensile tests for DC06 steel sheet samples with different alignment to the rolling direction. As yield function/plastic potential, simple orthotropic quadratic Hill (1948) or non-quadratic Karafillis-Boyce (1993) stress functions are utilized. According to the non-associated flow rule approach, the anisotropy parameters of the yield function are adjusted to the data related to the evolution of the yield stress ratios while the anisotropy parameters of the plastic potential are adjusted to the data related to evolution of the r-values. Considering these yield stress ratio and r-value evolutions, based on the principle of the plastic work equivalence, the anisotropy parameters of the utilized stress functions are set as 4^th order polynomial functions of the equivalent plastic strain. For the considered DC06 sheet sample, the analyzed non-associated formulations result in acceptable predictions of the yield stress and r-value directional dependences and their evolutions with on-going deformation. Besides, considering yield function/plastic potential contours, evolution of the plastic potential is more pronounced especially for the quadratic plastic potential function. Therefore, for the considered material, the presented results clearly indicate advantage of the non-associated over the associated flow rule approach. Furthermore, in the constitutive formulation, for the equivalent plastic strain as the internal variable that governs the evolution of the yield function/plastic potential, evolution equation consistent with the same principle of plastic work equivalence is utilized.

In the present paper, the algorithmic formulations of the analyzed non-associated constitutive descriptions with yield function/plastic potential evolution are developed based on the implicit Euler backward integration procedure. Application of this procedure results in the system of four non-linear algebraic equations (tensorial and scalar) that present incremental form of the constitutive equations. For the obtained system, the solution procedure based on the iterative Newton-Raphson procedure is developed. Considering data for DC06 sheet sample, accuracy of the derived stress integration procedures has been estimated by calculating iso-error maps for three representative stress points. For the analyzed non-associated formulations with yield function/plastic potential evolution, the error magnitudes are comparable to the error magnitudes obtained for the associated/non-associated formulations based on the same orthotropic stress functions and the isotropic hardening concept. It has been revealed that although there is a difference in maps configuration there is no considerable difference in error magnitudes for the formulations based on the different flow rule approach and the quadratic or the non-quadratic stress function. Furthermore, based on the obtained results it can be stated that for the analyzed formulations, a reduction of error magnitude can be expected when the strain size increment is reduced.

References

• An, Y.G., Vegter, H., Melzer, S., Triguero, P.R. (2013). Evolution of the plastic anisotropy with straining and its implication on formability for sheet metals, Journal of Materials Processing Technology 213: 1419-1425.
• Aretz, H. (2005). A non-quadratic plane stress yield function for orthotropic sheet metals, Journal of Materials Processing Technology 168: 1-9.
• Banabic D., Kuwabara T., Balan T., Comsa D. S., Julean D. (2003). Non-Quadratic yield criterion for orthotropic sheet metals under plane-stress conditions, The International Journal of Mechanical Sciences 45: 797-811.
• Banabic, D., Aretz, H., Comsa, D.S., Paraianu, L. (2005). An improved analytical description of orthotropy in metallic sheets, International Journal of Plasticity 21: 493-512.
• Barlat, F., Aretz, H., Yoon, J.W., Karabin, M.E., Brem, J.C., Dick, R.E. (2005). Linear transformation-based anisotropic yield functions, International Journal of Plasticity 21: 1009-1039.
• Barlat, F., Brem, J.C., Yoon, J.W., Chung, K., Dick, R.E., Lege, D.J., Pourboghrat, F., Choi, S.H., Chu, E. (2003). Plane stress yield function for aluminum alloy sheets - part 1: theory, International Journal of Plasticity 19: 1297-1319.
• Barlat, F., Lege, D.J., Brem, J.C. (1991). A six-component yield function for anisotropic metals, International Journal of Plasticity 7: 693-712.
• Barlat, F., Lian, J. (1989). Plastic behavior and stretchability of sheet metals. Part I: A yield function for orthotropic sheets under plane stress conditions, International Journal of Plasticity 5: 51-66.
• Barlat, F., Maeda, Y., Chung, K., Yanagawa, M., Brem, J.C., Hayashida, Y., Lege, D.J., Matsui, K., Murtha, S.J., Hattori, S., Becker, R.C., Makosey, S. (1997). Yield function development for aluminum alloy sheets, Journal of the Mechanics and Physics of Solids 45: 1727-1763.
• Cvitanić, V., Kovačić, M., Vladislavić, A. (2016). Numerical analysis of accuracy for evolutionary anisotropic plasticity models, Engineering Review 36: 255-267.
• Cvitanić, V., Vlak, F., Lozina, Ž. (2008). A finite element formulation based on non-associated plasticity for sheet metal forming, International Journal of Plasticity 24: 646-687.
• Hill, R. (1948). Theory of yielding and plastic flow of anisotropic metals, Proceedings of the Royal Society A 193: 281-297.
• Hill, R. (1979). Theoretical plasticity of textured aggregates, Mathematical Proceedings of the Cambridge Philosophical Society 55: 179-191.
• Hill, R. (1990). Constitutive modeling of orthotropic plasticity in sheet metals, Journal of the Mechanics and Physics of Solids 38: 405-417.
• Hill, R. (1993). A user-friendly theory of orthotropic plasticity in sheet metals, International Journal of Mechanical Sciences 35: 19-25.
• Karafillis, A.P., Boyce, M. (1993). A general anisotropic yield criterion using bounds and a transformation weighting tensor, Journal of the Mechanics and Physics of Solids 41: 1859-1886.
• Lademo, O.-G., Hopperstad, O.S., Langseth, M. (1999). An evaluation of yield criteria and flow rules for aluminium alloys, International Journal of Plasticity 15: 191-208.
• Mroz, Z. (1963). Non-associated flow laws in plasticity, J. de Mècanique 2: 21-42.
• Park, T., Chung, K. (2012). Non-associated flow rule with symmetric stiffness modulus for isotropic-kinematic hardening and its application for earing in circular cup drawing, International Journal of Solids and Structures 49: 3582-3593.
• Runesson, K., Mroz, Z. (1989). A note on non-associated plastic flow rules, International Journal of Plasticity 5: 639-658.
• Safaei, M., Lee, M.-G., Zang S.-I., De Waele, W. (2014). An evolutionary anisotropic model for sheet metals based on non-associated flow rule approach, Computational Materials Science 81: 15-29.
• Safaei, M., Zang, S.-I., Lee, M.-G., De Waele, W. (2013). Evaluation of anisotropic constitutive models: Mixed anisotropic hardening and non-associated flow rule approach, International Journal of Mechanical Sciences 73: 53-68.
• Simo, J.C., Hughes T.J.R. (1988). Elastoplasticity and Viscoplasticity - Computational Aspects, Springer-Verlag.
• Stoughton, T.B. (2002). A non-associated flow rule for sheet metal forming, International Journal of Plasticity 18: 687-714.
• Stoughton, T.B. Yoon, J.W., (2004). A pressure-sensitive yield criterion under a non-associated flow rule for sheet metal forming, International Journal of Plasticity 20: 705-731.
• Stoughton, T.B., Yoon, J.W. (2006). Review of Drucker’s postulate and the issue of plastic stability in metal forming, International Journal of Plasticity 22: 391-433.
• Taherizadeh, A., Green, D.E., Ghaei, A., Yoon, J.W. (2010). A non-associated constitutive model with mixed iso-kinematic hardening for finite element simulation of sheet metal forming, International Journal of Plasticity 26: 288-309.
• Taherizadeh, A., Green, D.E., Yoon, J.W. (2011). Evaluation of advanced anisotropic models with mixed hardening for general associated and non-associated flow metal plasticity, International Journal of Plasticity 27: 1781-1802.
• Vrh, M., Halilovič, M., Starman, B., Štok, B., Comsa, D.-S., Banabic, D. (2014). Capability of the BBC2008 yield criterion in predicting the earing profile in cup deep drawing simulations, European Journal of Mechanics; A/Solids 45: 59-74.
• Yoon, J.W., Barlat, F., Chung, K, Pourboghrat, F., Yang, D.Y. (2000). Earing prediction based on aysmmetric nonquadratic yield function, International Journal of Plasticity 216: 1075-1104.
• Yoon, J.W., Barlat, F., Chung, K, Pourboghrat, F., Yang, D.Y. (1998). Influence of initial back stress on the earing prediction of drawn cups for planar anisotropic aluminum sheets, Journal of Materials Processing Technology 80-81: 433-437.
• Yoon, J.W., Barlat, F., Dick, R.E., Chung, K., Kang, T.J. (2004). Plane stress yield function for aluminum alloy sheets - part II: FE formulation and its implementation, International Journal of Plasticity 20: 495-522.
• Yoon, J.W., Barlat, F., Dick, R.E., Karabin, M.E. (2006). Prediction of six or eight ears in a drawn cup based on a new anisotropic yield function, International Journal of Plasticity 22: 174-193.
• Yoon, J.W., Yang, D.Y., Chung, K. (1999). Elasto-plastic finite element method based on incremental deformation theory and continuum based shell elements for planar anisotropic sheet materials, Computer Methods in Applied Mechanics and Engineering 174: 23-56.
• Zamiri, A. and Pourboghrat, F. (2007). Characterization and development of an evolutionary yield function for the superconducting niobium sheet, International Journal of Solids and Structures 44: 8627-8647.

Publication Dates

• Publication in this collection
  Oct 2017

History

• Received
  13 Oct 2016
• Reviewed
  09 Mar 2017
• Accepted
  31 May 2017