- Email: [email protected]
Explicit incremental-update algorithm for modeling crystal elasto-viscoplastic response in finite element simulation
Available online at www.sciencedirect.com 8CICNGR
Science Press
d
DIRECT'
Transactions of Nonferrous Metals Society of China
Trans. Nonferrous Met. SOC.China 16(2006) s624-s630 www.csu.edu.cn/ysxb/
Explicit incremental-update algorithm for modeling crystal elasto-viscoplastic response in finite element simulation LI Hong-wei(%%fb), YANG He(#& $?),SUN Zhi-chao(.S/l&&B) College of Materials Science and Engineering, Northwestern Polytechnical University, Xi'an 7 10072, China Received 10 April 2006; accepted 25 April 2006 Abstract: Computational stability and efficiency are the key problems for numerical modeling of crystal plasticity, which will limit its development and application in finite element (FE) simulation evidently. Since implicit iterative algorithms are inefficient and have difficulty to determine initial values, an explicit incremental-update algorithm for the elasto-viscoplastic constitutive relation was developed in the intermediate frame by using the second Piola-Kirchoff (P-K) stress and Green stain. The increment of stress and slip resistance were solved by a calculation loop of linear equations sets. The reorientation of the crystal as well as the elastic strain can be obtained from a polar decomposition of the elastic deformation gradient. User material subroutine VUMAT was developed to combine crystal elasto-viscoplastic constitutive model with ABAQUSiExplicit. Numerical studies were performed on a cubic upset model with OFHC material (FCC crystal). The comparison of the numerical results with those obtained by implicit iterative algorithm and those From experiments demonstrates that the explicit algorithm is reliable. Furthermore, the effect rules of material anisotropy, rate sensitivity coefficient (RSC) and loading speeds on the deformation were studied. The numerical studies indicate that the explicit algorithm is suitable and efficient for large deformation analyses where anisotropy due to texture is important. Key words: crystal plasticity; explicit algorithm; constitutive relation; FE simulation; ABAQUSiExplicit
1 Introduction The microscopic constitutive relation (MCR) of crystal material is established on the base of crystal plasticity theory. It describes the mechanism of crystalline dislocation slip in plastic deformation of crystal. The combination of MCR with FEM, namely, crystal plastic finite element method (CPFEM), is an effective way to simulate and analyze crystal deformation in microscopic scales. CPFEM was proposed first by HILL and RICE[l], and the numerical method was established by PIERCE et a1[2]. Herefiom, a wide use of CPFEM has been achieved in the fields of material and mechanics. Since most materials are rate dependent, ratedependent crystal plasticity is widely used. Whereas, it is difficult to achieve the numerical application of rate-dependent crystal plasticity in FEM due to the high-order exponent hardening equations, which results
in the computational instability. To ensure a stable calculation, robust implicit iterative algorithms were introduced. KALIDINDI et a1[3] deduced a timeincremental crystal plasticity model and adopted Newton iteration method for the solution of the constitutive equations. To ensure computational convergence, additional limits were imposed on the stress increment. SARMA and ZACHARIA[4] provided an integration algorithm based on a time-incremental formula of the elastic deformation gradient. In the same way, the change of the elastic deformation gradient was limited to aid convergence. Three ordinary differential equations of the crystal plastic constitutive relation were derived by RASHID and NEMAT-NASSER[5], and a weighted residual method was adopted for the approximate solutions of these equations. RAPHANEL et a1[6] proposed an equilibrium search method based on the second-order RUNGE-KUTTA integration scheme, which achieved higher precision and a quadratic rate of convergence. However, it is still complicated. So, it can
Foundation item: Project(50335060) supported by the National Natural Science Foundation of China; Project (502255 18) supported by the National Science Fund for Distinguished Young Scholars of China; Project supported by the Scientific and Technological Innovation Foundation for Youth NPU Teachers Corresponding author: YANG He; Tel: +86-29-88495632; E-mail: [email protected]
LI Hong-wei, et a1 /Trans. Nonferrous Met. SOC.China 16(2006)
be concluded that plenty of iterations and the difficulty to determine initial iterative values make implicit algorithms inefficient and incompetent for simulations of large structure complex deformation process. Therefore, explicit algorithms are deduced to overcome these disadvantages. The algorithms are based on the linearization of the nonlinear constitutive relations. PIERCE et a1[7] proposed a forward-gradient method (i.e. tangent-modulus method), which is widely used in the subsequent numerical study of crystal plasticity. Whereas, it leads to a stiff system in calculation, which requires very small time increment in order to obtain the desired accuracy and stability. So, a plastic-predictor elasticcorrector method was proposed by NEMAT-NASSER et a1[8, 91, which solved the slip system shearing rates by the linearizationof the deformation rate equation, instead of rate-dependent hardening equation. So, the approach keeps the high-order exponent hardening equation from calculations. Nevertheless, it is necessary to determine whether the slip systems are active, and to distinguish three different regimes based on the number of the active slip systems for the solution of the slip system shearing rates, which makes the approach complicated and inefficient. MCGINTY and MCDOWELL[101 proposed a new approach to identify active slip systems prior to determining their shearing rate, based on which a semi-implicit integration scheme was proposed. It can be found that stability and efficiency of computation are the key problems for the numerical application of crystal plasticity in FEM. Except for algorithms, the frame in which the constitutive model is established does effects on computational efficiency as well. Most works are performed in current frame based on the constitutive relations in the form of objective stress rate. However, the vectors of slip systems in current frame are time-dependent, while, they are timeindependent in reference and intermediate frames. SO, DELANNAY et al[ll] established a model of MCR in reference frame, which made it convenient to realize numerical application of crystal plasticity. Herein, an explicit incremental-update model was proposed in the intermediate frame for modeling the elasto-viscoplastic constitutive response by using the second P-K stress and Green stain. Linear incremental equation of the stress increment was deduced by using Taylor series expansion of the rate-dependent hardening equation, The complete pivot Gaussian elimination method was adopted for the solution, which made the calculation simple and stable. The following conventions will be used in the paper. The fourth-order elastic tensor and the second-order identity tensor were represented by C and I respectively. Tensors or vectors will be denoted by boldfaced italic. The product and inner product of two tensors A and B
s625
will be denoted by AB and A *B,respectively. The vector product of two vectors a and b will be denoted by a@b. The double dot metrix of the forth-order tensor C and second-order tensor A will be denoted by CA. The product of tensor and vector v will be denoted by Av. Quantities of AB, Av, CA, u@b and A B have the fallowing component representation: ( A B k = A&, (Av)~ = Aikvk, (C:A)u= Cv&u, (a@b)u=aibj and (A B)u= A&.
2
Crystal kinematics and relations
constitutive
2.1 Crystal kinematics
The general kinematics of the elastic-plastic deformation of crystal at finite strains was given by HILL and RICE [l], and developed by HILL and HAVNER[121. HAVNER[131 provided a comprehensive account of the subject. The total deformation gradient of finite strain from the reference frame to the current frame, F,was defined by
F = -ax
(1)
ax
where X and x denote the reference and current particle positions, respectively. The kinematics of crystal finite deformation is hinted in Fig. 1.
"/
/
\
F BO
Fig. 1 Schematic diagram of single crystal kinematics
The hypothesis was adopted that material deformed plastically by crystalline slip, while the lattice underwent elastic deformation. So, the total deformation gradient has a multiplicative decomposition as
F =F * F P (2) where F P is the plastic deformation gradient due to shearing along crystallographic slip planes. The plastic
LI Hong-wei, et a1 /Trans. Nonferrous Met. SOC.China 16(2006)
s626
deformation by slip was assumed to occur at constant volume with the orientation of the lattice remaining unchanged, so, det F =1 (plastic incompressibility).F' is the elastic deformation gradient which includes any rigid rotation of the lattice, and can be obtained by F * = F F p - ' , detF*>O
(3)
The velocity gradient in current frame is L=kF-l
=k*F*-'+F*@PFP-'F*-' =L* +LP
(4)
The plastic velocity gradient zb relative to the intermediate frame is expressed as a linear combination of the slip system shearing rates, as follows:
r,
a
a
The dyadic product of the unit vector in the slip direction, g a l with the unit normal vector of the slip plane, hiia, gives the Schmid orientation tensor, F a , in the intermediate frame. Details of vectors of slip systems for FCC crystal with 12 { 111}<11O> type slip systems are presented in Table 1.
Vectors of slip system (1 1i ) [ i 101 (1 1i)[oi 11 (1 1i)[ioi] ( i 1qi To] ( i i1)[101] (i1i)[oi i] ( i i i)[i i o ] ( i i i)[oi 11 ( i ii ) [ i o i ] (1 11)[1lo] (1 i i)[ioi] (1 i i ) [ o i i ]
1 2 3 4 5
6 7 8 9 10 11
12
(7)
The corresponding stress measure is the second P-K stress, f,which is the elastic work conjugate to the Green strain, given in terms of the Cauchy stress T by *-T
=Yo
(9)
Sb4ZjF1
where rn is RSC, yo is the reference rate of shearing. The resolved shear stress is the component of the traction along the direction of slip, and is related to the second P-K stress and the Schmid tensor in the intermediate frame as
(
7" = F * T *T~
*) . (P
@I % a ) = (c*T*
). P
(10)
For metallic materials, the elastic stretch is usually infinitesimal. For this situation, the resolved shear stress Z n may, with impunity, be approximated by ra = T 8 . p a (1 1) The slip resistance s" is taken to evolve as
(13)
The hardening matrix, q @ , describes the latent hardening behavior of a crystallite, which is complicated as described in literature [14]. For the 12 {111}<110> slip systems of FCC crystals, q@ is given by
T*= c : E * (6) The strain measure is given by the Green strain tensor, E*, defined in terms of the elastic deformation gradient, F*, and the second-order identity tensor, I, as
T*= (det F * )F*-'TF
1'"
'Im
h@ = q @ P ' ( n o sum o n p )
2.2 Constitutive relations The elastic constitutive relations of crystalline material are presented in the intermediate frame:
E =- F F * - I 2
Za
where h@ is the rate of strain hardening on slip system a due to a shearing on the slip system p, and is related to a single slip hardening rate, h", and the hardening matrix, q@ ,as
Table 1 Details of slip systems for FCC crystal
Number
rates on the slip system a due to crystalline slip, and is given in terms of resolved shear stress z" and a slip resistance sa. Here, a rate-dependent slip model with the power law is employed.
(8)
The quantity,v, in Eqn.(S) is the plastic shearing
where q is the ratio of the latent hardening rate to the self-hardening rate given by ZHOU et al[ 151, and A is a 3x3 matrix hlly populated by ones. The model for the single slip hardening rate in Eqn.(l3) proposed by KALIDINDI[3] is adopted.
where ho, a and s, are the slip system hardening parameters.
3 Incremental-update algorithm Assumption is adopted that t denotes the time at the start of the step, and t = t +At denotes the time at the end
LI Hong-wei, et a1 /Trans. Nonferrous Met. SOC. China 16(2006)
of the step. An evolution equation for F as
is presented
(16) where the increment of shearing strain on a slip system is given by
Ay" =
(z)At
s627
C a = C : [(1/2)Ba]
(28)
Then, substituting Eqns.( 17) and (1 8) into Eqn.(26), we obtained
T f ( z ) = T t n- A t x
(17) Rewriting Eqn.(29) in incremental form as
Taylor series expansion of Eqn.(9), neglecting the high-order terms, is expressed as
According to Eqn.( 12), the incremental evolution of the slip system resistance is presented as
According to Eqn.( 1 l), we obtained ~ 7 ="AT* .
AT* = T *(T)- T*(f)= T * n - T * ( t ) - At x
(19)
Substitute Eqns.(3) and (7) into Eqn.(6),
As" = b t C h @ ( s ' ( ~ ) ) ~ B
Let
x = F p - T (zPT (z)F(z)Fp-I(z)
(21)
Then, substituting Eqn.( 16) into Eqn.(21), we obtained
X = F p - T (t)FT ( r ) F ( z ) F f l( t ) -
p-'(t)c
Fp-T ( f ) F T (z)F(z)F
AyaFaT -
Eqn.(3O) is solved for AT', keeping As" fixed at its best available estimate. Since all values at the time t are known, Eqn.(3O) is a linear equation set in terms ofAT' . So, the complete pivot Gaussian elimination method is adopted for the solution. Then, a new As" is obtained based on AT* and the old As" through Eqn.(31). Then the recalculation of AT* is performed on the base of the n e w h a . The loop is continued until both AT*and&" satisfy the criterions as follows: dT* =ATgk+i - A T * k , and dT ij < l o -4 Sg (k=O,l)...) dsa
Neglecting the 0
terms, Eqn.(22) may be
=hak+l
- h a k ,
I *I
and (dsa) <
max
(32) w3S,
(k=O,l,...; and a =1,2,...,12) where so is the initial value of the slip resistance.
(33)
rewritten as
F P - ~ ( ~ ) F T ( ~ ) F ( ~ ) F ~A--' (c ~A) =~ ~ B ~(23)
4 Procedures of algorithm
a
where
In the procedures, the time-independent slip systems mu ), as well as the deformation 'gradient, the plastic deformation gradient, the second P-K stress and the slip resistance at the start of the step, F ( t ),F ( t ) , T ' ( t )ands"(t), are known at first. We took all the initial values of s" (0)equal to so , and F (0) equal to the second-order identity tensor. Additionally, an estimate of the deformation gradient at time t,F ( z ) ,is given. The procedures are as follows: 1) Compute A, B" , T o mand C" using Eqns.(24), (25), (27) and (28), respectively; 2) Solve AT' and AP through the loop of ( Fa,
A = FP-T(t)FT(~)F(~)FP-'(t)
B" = AP"
+P u T A
(24) (25)
Finally, substituting Eqn.(23) into Eqn.(20), we obtained
T ' ( T ) =T" - C A y " C " where
T*' = C : [(1/2xA - I } ]
LI Hong-wei, et al /Trans. Nonferrous Met. SOC.China 16(2006)
s628
calculation by using Eqns.(30) and (31), and update T' (z) and su (z) ; 3) Update the plastic deformation gradient at the end of the step, F (z),by using Eqn.( 16). However, this does not ensure that det F (z) is equal to unity. So, it is necessary to normalize F (z) by
F;(z)=q(z)/@q
(34)
4) Compute the elastic deformation gradient at the end of the increment according to Eqn.(2), Fig. 2 Constraints situation in simple compression model
5) Compute the Cauchy stress in current frame by 1 F ( t )*T~(z) T(z)= det F * (z) 6) The polar decomposition ofF'(z) was used to compute the rotation, R' , of the crystal lattice during the deformation, F* =Rev* (37) Compute the reorientation of crystal lattice from +-T
sa(r)=R*Fa and m a ( z ) = R
-
ma
(38) 0
where sn (2) ,m a(z) represent a slip system in the current fiame at time t. It is convenient to implement the algorithm in FE simulation by developing user subroutine WMAT[ 161 on the platform of ABAQUS/Explicit.
5 Numerical results and analysis To ensure the subroutine is reliable prior to its wide use, the simulation of compression model with single element is adopted, as shown in Fig.2, with the material of the annealed oxygen free high conductivity (OFHC) copper, whose parameters are shown in Table 2. In the model of Fig.2, the length of side is 1 mm. The model is depressed with a speed of 20 mm.s-', and the total step time is 0.01 s. In addition, it is assumed that there is no friction on the faces. The responses of the axial stress and logarithmic axial strain obtained in the simulations are compared with those obtained by using the implicit iterative algorithm and those from experiment in the literature [3], as shown in Fig.3. The agreement is satisfying, which demonstrates the subroutine is reliable. Based on the above work, a cubic upset model with the side length of lOmm, meshed into 125 elements, was adopted for the simulations to do some research on the
0.2
0.4
0.6
0.8
Logarithmic axial strain,
Fig. 3 Comparison of axial stress-strain compression process
1 .O
1E I
responses in simple
effects of material anisotropy done on the deformation process. Nephograms of the distribution of stress in 3D view, as well as nephograms of the distribution of strain in y-z view at the different stages of the. deformation process are presented in Fig.4. The non-uniform distribution of the stress and strain, as well as the obvious drum shapes shown in Fig. 4 indicate that material anisotropy do much effect on the deformation, even though the load is uniaxial. In order to do some research on the effect rules of which loading speeds do on materials with different RSCs, the stress-strain responses of the material with the value 0.03 of RSC (m=0.03) under different loading speeds, and the stress-strain responses of materials with different RSCs under loading speed of 20 mm.s-' ( ~ 2 0 mm*s-') were investigated. The curves are shown in Fig.5 and Fig.6, respectively. In Fig.5, the loading speeds of 10, 50 and 100 m*ms-' are applied. It can be concluded that the axial stress increases with the increment of loading speed. Fig.6 shows that materials with RSCs of 0.02, 0.025 and
Table 2 Parameters of annealed OFHC copper Elastic modulus,
EIMPa 124 000
Poisson ratio 0.33
Reference value of slip rate, yoIs?
Rate sensitivity coefficient, m
Initial critical resolved shear stress, s"(O)/MPa
0.001
0.0 12
16.0
LI Hong-wei, et a1 /Trans. Nonferrous Met. SOC.China 16(2006)
s629
Fig. 4 Distribution of stress and strain with different deformation ratios 6: (a) Stress distribution, &lo%; (b) Strain distribution, 6=10%; (c) Stress distribution, 6=20%; (d) Strain distribution, 6=20%; (e) Stress distribution, 640%; (0 Strain distribution,640%
[
S
2
F
............
/.,
300
400
, n// ................ s ............... i ............... y
e
2
F-
300
b
9 2
.........
-.-m 2
-.-.-
I
0
200
* v)
I
I
I
I
0.2
0.4
0.6
0.8
.........
m=0.020
- m=0.025 --- m=0.030
100
I
I
I
I
1.0
Logarithmic axial strain, I E I Fig. 5 Stress-strain responses of material with m=0.03 under different loading speeds
0.03 undergo the deformation, and the axial stress increases with the increment of RSC. This rules agree well with the nature.
The simulations on the cubic upset model above were performed on the personal computer with the configuretion of Pentium 4, 2.60 GHz CPU and 2.59 GHz, 2.00 GB memory. Computational efficiency was
LI Hong-wei, et al /Trans. Nonferrous Met. SOC.China 16(2006)
s630
compared with the simulations by using the implicit algorithm in literature [ 171, which is shown in Table 3. Table 3 Comparison of computational efficiency between explicit and implicit algorithm Distance CPU time by CPU time by Number Of depressed using explicit using implicit element algorithmh algorithmh /mm 1 125
3.98 4.01
23 2 604
192 28 023
From Table 3, it can be seen that the explicit algorithm proposed is greatly efficient.
6 Conclusions The explicit incremental-update algorithm of crystal elasto-viscoplastic constitutive relation was presented by using the second P-K stress and Green stain in the intermediate frame. The increments of stress and the slip resistance are solved by a loop calculation of linear equations set, which is simplified on the base of the linearization of the nonlinear equations. User material subroutine WMAT is developed in the light of the algorithm to implement crystal elasto-viscoplastic constitutive model in FE simulation. The comparison of the numerical results with those obtained by using implicit iterative algorithm and those from experiment demonstrates the explicit algorithm in the paper is reliable. Further studies show that material anisotropy, RSC and loading speed do effects on the deformation, the axial stress increases with the increment of loading speed and RSC. The studies indicate that the explicit algorithm is reliable and efficient.
References [I]
HILL R, RICE J R. Constitutive analysis of elastic-plastic crystals at
arbitrary strain [J]. J Mech Phys Solids, 1972, 20(6): 401-413. PIERCE D, ASARO R J, NEEDLEEMAN A. Material rate dependence and localized deformation in crystalline solids [J]. Acta Metall, 1983,31(12): 1951-1976. KALIDINDI S R, BRONKHORST C A, ANAND L. Crystallographic texture evolution in bulk deformation process of fcc metals [J]. J Mech Phys Solids, 1992,40(3): 537-569. SA R M A G, ZACHARIA T. Integration algorithm for modeling the elasto-viscoplastic response of polycrystalline materials [J]. J Mech Phys Solids, 1999,47(6): 1219-1238. RASHID M M, NEMAT-NASSER S. A constitutive algorithm for rate- dependent crystal plasticity [J]. Computer Methods in Applied Mechanics and Engineering, 1992,94(2): 201- 228. RAF'HANEL J L, RAVICHANDRAN G, LEROY Y M. Threedimensional ratedependent crystal plasticity based on Runge-Kutta algorithms for update and consistent linearization [J]. International Journal of Solids and Structures, 2004,41(22-23): 5995-6021. PEIRCE D, SHIH C F, NEEDLEMAN A. A tangent modulus method for rate dependent solids [J]. Computers and Structures, 1984, 18(5): 875-887. NEMAT-NASSER S, OKINAKA T. A new computational approach to crystal plasticity: fcc single crystal [J]. Mechanics of Materials, 1996,24(1): 43-57. NEMAT-NASSER S, NI L, OKMAKA T. A constitutive model for fcc crystals with application to polycrystalline OFHC copper [J]. Mechanics of Materials, 1998, 30: 325-341. MCGINTY R D, MCDOWELL D L. A semi-implicit integration for rate independent finite crystal plasticity [J]. Int J Plasticity, 2006, 22(6): 996-1025. DELANNAY L, JACQUES P J, KALIDINDI. S R. Finite element modeling of crystal plasticity with grains shaped as truncated octahedrons [J]. Int J Plasticity, 2006. (in press) HILL R, HAVNER KS. Perspectives in the mechanics of elastoplastic crystal [J]. J Mech Phys Solids, 1982,30( 1-2): 5-22. HAVNER K S. Finite plastic deformation of crystalline solids [MI. New York: Cambridge, 1992. KUMAR A V, YANG C. Study of work hardening models for single crystals using three dimensional finite element analysis [J]. Int J Plasticity, 1999, 15(7): 737-754. ZHOU Y, NEALE K W, T6TH L S. A modified model for simulating latent hardening during the plastic deformation of ratedependent fcc polycrystals [J]. Int J Plasticity, 1993, 9(8): 96 1-978. ABAQUS User's Manual, 2004. Version 6.5, ABAQUS, tnc. LI H W, YANG H, SUN Z C. Research on the implicit algorithm and numeralization of microscopic constitutive relation [J]. Mater Sci Technol, 2006. (Accepted) (Edited by LI Xieng-qun)
Science Press
d
DIRECT'
Transactions of Nonferrous Metals Society of China
Trans. Nonferrous Met. SOC.China 16(2006) s624-s630 www.csu.edu.cn/ysxb/
Explicit incremental-update algorithm for modeling crystal elasto-viscoplastic response in finite element simulation LI Hong-wei(%%fb), YANG He(#& $?),SUN Zhi-chao(.S/l&&B) College of Materials Science and Engineering, Northwestern Polytechnical University, Xi'an 7 10072, China Received 10 April 2006; accepted 25 April 2006 Abstract: Computational stability and efficiency are the key problems for numerical modeling of crystal plasticity, which will limit its development and application in finite element (FE) simulation evidently. Since implicit iterative algorithms are inefficient and have difficulty to determine initial values, an explicit incremental-update algorithm for the elasto-viscoplastic constitutive relation was developed in the intermediate frame by using the second Piola-Kirchoff (P-K) stress and Green stain. The increment of stress and slip resistance were solved by a calculation loop of linear equations sets. The reorientation of the crystal as well as the elastic strain can be obtained from a polar decomposition of the elastic deformation gradient. User material subroutine VUMAT was developed to combine crystal elasto-viscoplastic constitutive model with ABAQUSiExplicit. Numerical studies were performed on a cubic upset model with OFHC material (FCC crystal). The comparison of the numerical results with those obtained by implicit iterative algorithm and those From experiments demonstrates that the explicit algorithm is reliable. Furthermore, the effect rules of material anisotropy, rate sensitivity coefficient (RSC) and loading speeds on the deformation were studied. The numerical studies indicate that the explicit algorithm is suitable and efficient for large deformation analyses where anisotropy due to texture is important. Key words: crystal plasticity; explicit algorithm; constitutive relation; FE simulation; ABAQUSiExplicit
1 Introduction The microscopic constitutive relation (MCR) of crystal material is established on the base of crystal plasticity theory. It describes the mechanism of crystalline dislocation slip in plastic deformation of crystal. The combination of MCR with FEM, namely, crystal plastic finite element method (CPFEM), is an effective way to simulate and analyze crystal deformation in microscopic scales. CPFEM was proposed first by HILL and RICE[l], and the numerical method was established by PIERCE et a1[2]. Herefiom, a wide use of CPFEM has been achieved in the fields of material and mechanics. Since most materials are rate dependent, ratedependent crystal plasticity is widely used. Whereas, it is difficult to achieve the numerical application of rate-dependent crystal plasticity in FEM due to the high-order exponent hardening equations, which results
in the computational instability. To ensure a stable calculation, robust implicit iterative algorithms were introduced. KALIDINDI et a1[3] deduced a timeincremental crystal plasticity model and adopted Newton iteration method for the solution of the constitutive equations. To ensure computational convergence, additional limits were imposed on the stress increment. SARMA and ZACHARIA[4] provided an integration algorithm based on a time-incremental formula of the elastic deformation gradient. In the same way, the change of the elastic deformation gradient was limited to aid convergence. Three ordinary differential equations of the crystal plastic constitutive relation were derived by RASHID and NEMAT-NASSER[5], and a weighted residual method was adopted for the approximate solutions of these equations. RAPHANEL et a1[6] proposed an equilibrium search method based on the second-order RUNGE-KUTTA integration scheme, which achieved higher precision and a quadratic rate of convergence. However, it is still complicated. So, it can
Foundation item: Project(50335060) supported by the National Natural Science Foundation of China; Project (502255 18) supported by the National Science Fund for Distinguished Young Scholars of China; Project supported by the Scientific and Technological Innovation Foundation for Youth NPU Teachers Corresponding author: YANG He; Tel: +86-29-88495632; E-mail: [email protected]
LI Hong-wei, et a1 /Trans. Nonferrous Met. SOC.China 16(2006)
be concluded that plenty of iterations and the difficulty to determine initial iterative values make implicit algorithms inefficient and incompetent for simulations of large structure complex deformation process. Therefore, explicit algorithms are deduced to overcome these disadvantages. The algorithms are based on the linearization of the nonlinear constitutive relations. PIERCE et a1[7] proposed a forward-gradient method (i.e. tangent-modulus method), which is widely used in the subsequent numerical study of crystal plasticity. Whereas, it leads to a stiff system in calculation, which requires very small time increment in order to obtain the desired accuracy and stability. So, a plastic-predictor elasticcorrector method was proposed by NEMAT-NASSER et a1[8, 91, which solved the slip system shearing rates by the linearizationof the deformation rate equation, instead of rate-dependent hardening equation. So, the approach keeps the high-order exponent hardening equation from calculations. Nevertheless, it is necessary to determine whether the slip systems are active, and to distinguish three different regimes based on the number of the active slip systems for the solution of the slip system shearing rates, which makes the approach complicated and inefficient. MCGINTY and MCDOWELL[101 proposed a new approach to identify active slip systems prior to determining their shearing rate, based on which a semi-implicit integration scheme was proposed. It can be found that stability and efficiency of computation are the key problems for the numerical application of crystal plasticity in FEM. Except for algorithms, the frame in which the constitutive model is established does effects on computational efficiency as well. Most works are performed in current frame based on the constitutive relations in the form of objective stress rate. However, the vectors of slip systems in current frame are time-dependent, while, they are timeindependent in reference and intermediate frames. SO, DELANNAY et al[ll] established a model of MCR in reference frame, which made it convenient to realize numerical application of crystal plasticity. Herein, an explicit incremental-update model was proposed in the intermediate frame for modeling the elasto-viscoplastic constitutive response by using the second P-K stress and Green stain. Linear incremental equation of the stress increment was deduced by using Taylor series expansion of the rate-dependent hardening equation, The complete pivot Gaussian elimination method was adopted for the solution, which made the calculation simple and stable. The following conventions will be used in the paper. The fourth-order elastic tensor and the second-order identity tensor were represented by C and I respectively. Tensors or vectors will be denoted by boldfaced italic. The product and inner product of two tensors A and B
s625
will be denoted by AB and A *B,respectively. The vector product of two vectors a and b will be denoted by a@b. The double dot metrix of the forth-order tensor C and second-order tensor A will be denoted by CA. The product of tensor and vector v will be denoted by Av. Quantities of AB, Av, CA, u@b and A B have the fallowing component representation: ( A B k = A&, (Av)~ = Aikvk, (C:A)u= Cv&u, (a@b)u=aibj and (A B)u= A&.
2
Crystal kinematics and relations
constitutive
2.1 Crystal kinematics
The general kinematics of the elastic-plastic deformation of crystal at finite strains was given by HILL and RICE [l], and developed by HILL and HAVNER[121. HAVNER[131 provided a comprehensive account of the subject. The total deformation gradient of finite strain from the reference frame to the current frame, F,was defined by
F = -ax
(1)
ax
where X and x denote the reference and current particle positions, respectively. The kinematics of crystal finite deformation is hinted in Fig. 1.
"/
/
\
F BO
Fig. 1 Schematic diagram of single crystal kinematics
The hypothesis was adopted that material deformed plastically by crystalline slip, while the lattice underwent elastic deformation. So, the total deformation gradient has a multiplicative decomposition as
F =F * F P (2) where F P is the plastic deformation gradient due to shearing along crystallographic slip planes. The plastic
LI Hong-wei, et a1 /Trans. Nonferrous Met. SOC.China 16(2006)
s626
deformation by slip was assumed to occur at constant volume with the orientation of the lattice remaining unchanged, so, det F =1 (plastic incompressibility).F' is the elastic deformation gradient which includes any rigid rotation of the lattice, and can be obtained by F * = F F p - ' , detF*>O
(3)
The velocity gradient in current frame is L=kF-l
=k*F*-'+F*@PFP-'F*-' =L* +LP
(4)
The plastic velocity gradient zb relative to the intermediate frame is expressed as a linear combination of the slip system shearing rates, as follows:
r,
a
a
The dyadic product of the unit vector in the slip direction, g a l with the unit normal vector of the slip plane, hiia, gives the Schmid orientation tensor, F a , in the intermediate frame. Details of vectors of slip systems for FCC crystal with 12 { 111}<11O> type slip systems are presented in Table 1.
Vectors of slip system (1 1i ) [ i 101 (1 1i)[oi 11 (1 1i)[ioi] ( i 1qi To] ( i i1)[101] (i1i)[oi i] ( i i i)[i i o ] ( i i i)[oi 11 ( i ii ) [ i o i ] (1 11)[1lo] (1 i i)[ioi] (1 i i ) [ o i i ]
1 2 3 4 5
6 7 8 9 10 11
12
(7)
The corresponding stress measure is the second P-K stress, f,which is the elastic work conjugate to the Green strain, given in terms of the Cauchy stress T by *-T
=Yo
(9)
Sb4ZjF1
where rn is RSC, yo is the reference rate of shearing. The resolved shear stress is the component of the traction along the direction of slip, and is related to the second P-K stress and the Schmid tensor in the intermediate frame as
(
7" = F * T *T~
*) . (P
@I % a ) = (c*T*
). P
(10)
For metallic materials, the elastic stretch is usually infinitesimal. For this situation, the resolved shear stress Z n may, with impunity, be approximated by ra = T 8 . p a (1 1) The slip resistance s" is taken to evolve as
(13)
The hardening matrix, q @ , describes the latent hardening behavior of a crystallite, which is complicated as described in literature [14]. For the 12 {111}<110> slip systems of FCC crystals, q@ is given by
T*= c : E * (6) The strain measure is given by the Green strain tensor, E*, defined in terms of the elastic deformation gradient, F*, and the second-order identity tensor, I, as
T*= (det F * )F*-'TF
1'"
'Im
h@ = q @ P ' ( n o sum o n p )
2.2 Constitutive relations The elastic constitutive relations of crystalline material are presented in the intermediate frame:
E =- F F * - I 2
Za
where h@ is the rate of strain hardening on slip system a due to a shearing on the slip system p, and is related to a single slip hardening rate, h", and the hardening matrix, q@ ,as
Table 1 Details of slip systems for FCC crystal
Number
rates on the slip system a due to crystalline slip, and is given in terms of resolved shear stress z" and a slip resistance sa. Here, a rate-dependent slip model with the power law is employed.
(8)
The quantity,v, in Eqn.(S) is the plastic shearing
where q is the ratio of the latent hardening rate to the self-hardening rate given by ZHOU et al[ 151, and A is a 3x3 matrix hlly populated by ones. The model for the single slip hardening rate in Eqn.(l3) proposed by KALIDINDI[3] is adopted.
where ho, a and s, are the slip system hardening parameters.
3 Incremental-update algorithm Assumption is adopted that t denotes the time at the start of the step, and t = t +At denotes the time at the end
LI Hong-wei, et a1 /Trans. Nonferrous Met. SOC. China 16(2006)
of the step. An evolution equation for F as
is presented
(16) where the increment of shearing strain on a slip system is given by
Ay" =
(z)At
s627
C a = C : [(1/2)Ba]
(28)
Then, substituting Eqns.( 17) and (1 8) into Eqn.(26), we obtained
T f ( z ) = T t n- A t x
(17) Rewriting Eqn.(29) in incremental form as
Taylor series expansion of Eqn.(9), neglecting the high-order terms, is expressed as
According to Eqn.( 12), the incremental evolution of the slip system resistance is presented as
According to Eqn.( 1 l), we obtained ~ 7 ="AT* .
AT* = T *(T)- T*(f)= T * n - T * ( t ) - At x
(19)
Substitute Eqns.(3) and (7) into Eqn.(6),
As" = b t C h @ ( s ' ( ~ ) ) ~ B
Let
x = F p - T (zPT (z)F(z)Fp-I(z)
(21)
Then, substituting Eqn.( 16) into Eqn.(21), we obtained
X = F p - T (t)FT ( r ) F ( z ) F f l( t ) -
p-'(t)c
Fp-T ( f ) F T (z)F(z)F
AyaFaT -
Eqn.(3O) is solved for AT', keeping As" fixed at its best available estimate. Since all values at the time t are known, Eqn.(3O) is a linear equation set in terms ofAT' . So, the complete pivot Gaussian elimination method is adopted for the solution. Then, a new As" is obtained based on AT* and the old As" through Eqn.(31). Then the recalculation of AT* is performed on the base of the n e w h a . The loop is continued until both AT*and&" satisfy the criterions as follows: dT* =ATgk+i - A T * k , and dT ij < l o -4 Sg (k=O,l)...) dsa
Neglecting the 0
terms, Eqn.(22) may be
=hak+l
- h a k ,
I *I
and (dsa) <
max
(32) w3S,
(k=O,l,...; and a =1,2,...,12) where so is the initial value of the slip resistance.
(33)
rewritten as
F P - ~ ( ~ ) F T ( ~ ) F ( ~ ) F ~A--' (c ~A) =~ ~ B ~(23)
4 Procedures of algorithm
a
where
In the procedures, the time-independent slip systems mu ), as well as the deformation 'gradient, the plastic deformation gradient, the second P-K stress and the slip resistance at the start of the step, F ( t ),F ( t ) , T ' ( t )ands"(t), are known at first. We took all the initial values of s" (0)equal to so , and F (0) equal to the second-order identity tensor. Additionally, an estimate of the deformation gradient at time t,F ( z ) ,is given. The procedures are as follows: 1) Compute A, B" , T o mand C" using Eqns.(24), (25), (27) and (28), respectively; 2) Solve AT' and AP through the loop of ( Fa,
A = FP-T(t)FT(~)F(~)FP-'(t)
B" = AP"
+P u T A
(24) (25)
Finally, substituting Eqn.(23) into Eqn.(20), we obtained
T ' ( T ) =T" - C A y " C " where
T*' = C : [(1/2xA - I } ]
LI Hong-wei, et al /Trans. Nonferrous Met. SOC.China 16(2006)
s628
calculation by using Eqns.(30) and (31), and update T' (z) and su (z) ; 3) Update the plastic deformation gradient at the end of the step, F (z),by using Eqn.( 16). However, this does not ensure that det F (z) is equal to unity. So, it is necessary to normalize F (z) by
F;(z)=q(z)/@q
(34)
4) Compute the elastic deformation gradient at the end of the increment according to Eqn.(2), Fig. 2 Constraints situation in simple compression model
5) Compute the Cauchy stress in current frame by 1 F ( t )*T~(z) T(z)= det F * (z) 6) The polar decomposition ofF'(z) was used to compute the rotation, R' , of the crystal lattice during the deformation, F* =Rev* (37) Compute the reorientation of crystal lattice from +-T
sa(r)=R*Fa and m a ( z ) = R
-
ma
(38) 0
where sn (2) ,m a(z) represent a slip system in the current fiame at time t. It is convenient to implement the algorithm in FE simulation by developing user subroutine WMAT[ 161 on the platform of ABAQUS/Explicit.
5 Numerical results and analysis To ensure the subroutine is reliable prior to its wide use, the simulation of compression model with single element is adopted, as shown in Fig.2, with the material of the annealed oxygen free high conductivity (OFHC) copper, whose parameters are shown in Table 2. In the model of Fig.2, the length of side is 1 mm. The model is depressed with a speed of 20 mm.s-', and the total step time is 0.01 s. In addition, it is assumed that there is no friction on the faces. The responses of the axial stress and logarithmic axial strain obtained in the simulations are compared with those obtained by using the implicit iterative algorithm and those from experiment in the literature [3], as shown in Fig.3. The agreement is satisfying, which demonstrates the subroutine is reliable. Based on the above work, a cubic upset model with the side length of lOmm, meshed into 125 elements, was adopted for the simulations to do some research on the
0.2
0.4
0.6
0.8
Logarithmic axial strain,
Fig. 3 Comparison of axial stress-strain compression process
1 .O
1E I
responses in simple
effects of material anisotropy done on the deformation process. Nephograms of the distribution of stress in 3D view, as well as nephograms of the distribution of strain in y-z view at the different stages of the. deformation process are presented in Fig.4. The non-uniform distribution of the stress and strain, as well as the obvious drum shapes shown in Fig. 4 indicate that material anisotropy do much effect on the deformation, even though the load is uniaxial. In order to do some research on the effect rules of which loading speeds do on materials with different RSCs, the stress-strain responses of the material with the value 0.03 of RSC (m=0.03) under different loading speeds, and the stress-strain responses of materials with different RSCs under loading speed of 20 mm.s-' ( ~ 2 0 mm*s-') were investigated. The curves are shown in Fig.5 and Fig.6, respectively. In Fig.5, the loading speeds of 10, 50 and 100 m*ms-' are applied. It can be concluded that the axial stress increases with the increment of loading speed. Fig.6 shows that materials with RSCs of 0.02, 0.025 and
Table 2 Parameters of annealed OFHC copper Elastic modulus,
EIMPa 124 000
Poisson ratio 0.33
Reference value of slip rate, yoIs?
Rate sensitivity coefficient, m
Initial critical resolved shear stress, s"(O)/MPa
0.001
0.0 12
16.0
LI Hong-wei, et a1 /Trans. Nonferrous Met. SOC.China 16(2006)
s629
Fig. 4 Distribution of stress and strain with different deformation ratios 6: (a) Stress distribution, &lo%; (b) Strain distribution, 6=10%; (c) Stress distribution, 6=20%; (d) Strain distribution, 6=20%; (e) Stress distribution, 640%; (0 Strain distribution,640%
[
S
2
F
............
/.,
300
400
, n// ................ s ............... i ............... y
e
2
F-
300
b
9 2
.........
-.-m 2
-.-.-
I
0
200
* v)
I
I
I
I
0.2
0.4
0.6
0.8
.........
m=0.020
- m=0.025 --- m=0.030
100
I
I
I
I
1.0
Logarithmic axial strain, I E I Fig. 5 Stress-strain responses of material with m=0.03 under different loading speeds
0.03 undergo the deformation, and the axial stress increases with the increment of RSC. This rules agree well with the nature.
The simulations on the cubic upset model above were performed on the personal computer with the configuretion of Pentium 4, 2.60 GHz CPU and 2.59 GHz, 2.00 GB memory. Computational efficiency was
LI Hong-wei, et al /Trans. Nonferrous Met. SOC.China 16(2006)
s630
compared with the simulations by using the implicit algorithm in literature [ 171, which is shown in Table 3. Table 3 Comparison of computational efficiency between explicit and implicit algorithm Distance CPU time by CPU time by Number Of depressed using explicit using implicit element algorithmh algorithmh /mm 1 125
3.98 4.01
23 2 604
192 28 023
From Table 3, it can be seen that the explicit algorithm proposed is greatly efficient.
6 Conclusions The explicit incremental-update algorithm of crystal elasto-viscoplastic constitutive relation was presented by using the second P-K stress and Green stain in the intermediate frame. The increments of stress and the slip resistance are solved by a loop calculation of linear equations set, which is simplified on the base of the linearization of the nonlinear equations. User material subroutine WMAT is developed in the light of the algorithm to implement crystal elasto-viscoplastic constitutive model in FE simulation. The comparison of the numerical results with those obtained by using implicit iterative algorithm and those from experiment demonstrates the explicit algorithm in the paper is reliable. Further studies show that material anisotropy, RSC and loading speed do effects on the deformation, the axial stress increases with the increment of loading speed and RSC. The studies indicate that the explicit algorithm is reliable and efficient.
References [I]
HILL R, RICE J R. Constitutive analysis of elastic-plastic crystals at
arbitrary strain [J]. J Mech Phys Solids, 1972, 20(6): 401-413. PIERCE D, ASARO R J, NEEDLEEMAN A. Material rate dependence and localized deformation in crystalline solids [J]. Acta Metall, 1983,31(12): 1951-1976. KALIDINDI S R, BRONKHORST C A, ANAND L. Crystallographic texture evolution in bulk deformation process of fcc metals [J]. J Mech Phys Solids, 1992,40(3): 537-569. SA R M A G, ZACHARIA T. Integration algorithm for modeling the elasto-viscoplastic response of polycrystalline materials [J]. J Mech Phys Solids, 1999,47(6): 1219-1238. RASHID M M, NEMAT-NASSER S. A constitutive algorithm for rate- dependent crystal plasticity [J]. Computer Methods in Applied Mechanics and Engineering, 1992,94(2): 201- 228. RAF'HANEL J L, RAVICHANDRAN G, LEROY Y M. Threedimensional ratedependent crystal plasticity based on Runge-Kutta algorithms for update and consistent linearization [J]. International Journal of Solids and Structures, 2004,41(22-23): 5995-6021. PEIRCE D, SHIH C F, NEEDLEMAN A. A tangent modulus method for rate dependent solids [J]. Computers and Structures, 1984, 18(5): 875-887. NEMAT-NASSER S, OKINAKA T. A new computational approach to crystal plasticity: fcc single crystal [J]. Mechanics of Materials, 1996,24(1): 43-57. NEMAT-NASSER S, NI L, OKMAKA T. A constitutive model for fcc crystals with application to polycrystalline OFHC copper [J]. Mechanics of Materials, 1998, 30: 325-341. MCGINTY R D, MCDOWELL D L. A semi-implicit integration for rate independent finite crystal plasticity [J]. Int J Plasticity, 2006, 22(6): 996-1025. DELANNAY L, JACQUES P J, KALIDINDI. S R. Finite element modeling of crystal plasticity with grains shaped as truncated octahedrons [J]. Int J Plasticity, 2006. (in press) HILL R, HAVNER KS. Perspectives in the mechanics of elastoplastic crystal [J]. J Mech Phys Solids, 1982,30( 1-2): 5-22. HAVNER K S. Finite plastic deformation of crystalline solids [MI. New York: Cambridge, 1992. KUMAR A V, YANG C. Study of work hardening models for single crystals using three dimensional finite element analysis [J]. Int J Plasticity, 1999, 15(7): 737-754. ZHOU Y, NEALE K W, T6TH L S. A modified model for simulating latent hardening during the plastic deformation of ratedependent fcc polycrystals [J]. Int J Plasticity, 1993, 9(8): 96 1-978. ABAQUS User's Manual, 2004. Version 6.5, ABAQUS, tnc. LI H W, YANG H, SUN Z C. Research on the implicit algorithm and numeralization of microscopic constitutive relation [J]. Mater Sci Technol, 2006. (Accepted) (Edited by LI Xieng-qun)