- Email: [email protected]
Finite element simulation of large-strain single-crystal viscoplasticity: An investigation of various hardening relations
Computational Materials Science 81 (2014) 386–396
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Finite element simulation of large-strain single-crystal viscoplasticity: An investigation of various hardening relations T.M. Povall ⇑, A.T. McBride, B.D. Reddy Centre for Research in Computational and Applied Mechanics, University of Cape Town, Private Bag, 7701 Rondebosch, South Africa
a r t i c l e
i n f o
Article history: Received 6 May 2013 Received in revised form 31 July 2013 Accepted 21 August 2013 Available online 18 September 2013 Keywords: Single crystals Large strain Plasticity Viscoplasticity Hardening Finite element method
a b s t r a c t Hardening relations describe the increase in resistance to deformation during plastic flow. Three hardening relations are compared here in the context of conventional large-strain single-crystal viscoplasticity. The first is an isotropic hardening relation. The second is a hardening relation that is expressed as an ordinary differential equation in the slip resistance. The third is a new relation, originally developed in the context of gradient crystal plasticity, in which the slip resistance is expressed explicitly in terms of the accumulated slip on each slip system. The numerical solution of the governing equations is found using the finite element method coupled with a predictor–corrector type algorithm. The features of the hardening relations are elucidated using a series of numerical benchmark problems. The parameters for the hardening relations are calibrated using a model problem. Various crystal structures are investigated, including single- and double slip, and face-centred cubic crystals. The hardening relations are compared and their relative features discussed. Ó 2013 Elsevier B.V. All rights reserved.
1. Introduction The description of plastic behaviour is characterised by (a) a yield criterion; (b) a hardening relation; and (c) a flow rule for plastic strain rates or equivalent quantities. A variety of hardening relations have been proposed for singlecrystal plasticity. Popular models include a hardening relation that is expressed as an ordinary differential equation in the flow resistance [15,9]. A new hardening relation has recently been proposed by Gurtin and Reddy [8] in the context of strain-gradient single-crystal plasticity. This relation takes into account the interactions between slip systems. The purpose of this paper is to compare the fully-interactive hardening relation [8] with an implicit relation and an isotropic relation. The fully-interactive relation was originally proposed in the context of small-strain rate-independent single-crystal gradient plasticity. The relation is extended here to the conventional largestrain viscoplastic setting. The isotropic relation is that defined by Steinmann and Stein [21] and is dependent on the sum of the accumulated slips, while the implicit relation used is that given in [9]. The finite element method (FEM) is used to solve the system of governing equations. The FEM implementation for conventional single-crystal plasticity is a well-established field with contributions from Schröder and Miehe [18], Miehe et al. [13], Steinmann and Stein [21], and Anand and Kothari [1], to mention some impor⇑ Corresponding author. Tel.: +27 (0)21 650 3817. E-mail addresses: [email protected] (T.M. Povall), [email protected] (A.T. McBride), [email protected] (B.D. Reddy). 0927-0256/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.commatsci.2013.08.043
tant examples. A predictor–corrector algorithm, similar to that described by Steinmann and Stein [21], is used in conjunction with a consistent algorithmic tangent to accurately and efficiently describe the plastic evolution. Two model problems are used to compare the hardening relations. The first involves the indentation of an elastoplastic material by a rigid spherical indenter. The second problem is the shear of a bar subjected to cyclic deformations. The various hardening relations are examined in the indentation problem for three types of crystal structures: a single slip system, a double slip system and a face-centred cubic (FCC) crystal. The shear problem examines the hardening relations using a double slip-system crystal and a FCC crystal. This paper is arranged as follows. Section 2 gives an overview of the notation, kinematics and the governing equations. Section 3 provides a short derivation of the reduced dissipation inequality. A comparison of rate-independent and rate-dependent formulations is given in Section 4. Section 5 details the three hardening relations. An explanation of the calibration process is given in Section 6.1. The results of the numerical simulation of the indentation and shear problems are discussed in Sections 6.2 and 6.3, respectively. A summary of the results and the conclusion are given in Section 7. 2. Kinematics and governing equations 2.1. Kinematics A motion of a body B is a smooth function v that assigns to each material point X at time t a point in Euclidean space, x = v(X, t). The point x is referred to as the spatial (or current) point occupied by X
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
at time t. The body B in the reference configuration at time t0 is denoted by B0. The displacement u is defined by
uðX; tÞ ¼ vðX; tÞ X:
ð1Þ
The velocity v is defined as the time rate of change of the displacement, i.e.
_ v ðX; tÞ ¼ u:
ð2Þ
The deformation gradient F is defined by
F :¼ $v ¼
@ vi ei EA ; @X A
ð3Þ
where the Jacobian J = det F is positive and ei and EA represent bases in the current and reference configurations, respectively. The Kröner–Lee decomposition multiplicatively splits the deformation gradient into an elastic and a plastic part ([10,11]), denoted Fe and Fp respectively, and is given by
F ¼ Fe Fp :
det Fe > 0 and
det Fp ¼ 1:
ð5Þ
The multiplicative split motivates the introduction of an intermediate configuration. From the elastic deformation gradient, the right elastic Cauchy–Green tensor is defined by
Ce ¼ FeT Fe :
ð6Þ
The spatial velocity gradient L is the tensor field defined by
@v i 1 ¼ ei ej ¼ Le þ Fe Lp Fe ; @xj
ð7Þ
1
1
and Lp ¼ F_ p Fp :
ð8Þ
Note that in (7) the plastic distortion rate Lp is pushed forward from the intermediate (lattice) configuration into the current configuration using Fe. 2.2. Slip systems and the single-crystal hypothesis
L ¼
X
divr þ b ¼ 0;
ð11Þ
where the Cauchy stress r is a symmetric tensor and b is the body force per unit current volume. The divergence of a second order tensor A is defined as divA = @Aij/@xjei. The resolved shear stress sa is an important measure in crystal plasticity. It is defined as the Cauchy stress resolved on slip plane a in the current configuration:
a; sa ¼ sa rm
ð12Þ
a are the slip directions and normals pushed forward where sa and m into the spatial configuration from the lattice configuration, i.e.
a ¼ Fe T ma : and m
sa ¼ Fe sa
ð13Þ
3. Free energy and dissipation inequality The dissipation inequality arises from the second law of thermodynamics [6] and is given by
X 1 J 1 Sp : C_ e þ sa ma J1 w_ P 0; 2 a
where w is the free energy and Sp is a stress measure similar to the second Piola–Kirchhoff stress, but defined in the intermediate configuration as
a
m s m :
ð15Þ
It is assumed that the free energy is dependent on both the elastic and the plastic deformation. Furthermore, it is assumed that any hardening mechanism is uncoupled from the elastic deformation. The free energy function is thus of the form
w ¼ we ðCe Þ þ wp ðna Þ;
ð16Þ
where na is a set of internal scalar plastic variables. The time derivative of the free energy is thus
@we ðCe Þ _ e X@wp ðna Þ _ a w_ ¼ n : :C þ @na @Ce a |fflfflfflffl{zfflfflfflffl}
ð17Þ
The dissipation inequality becomes, on substitution of (18),
J 1
X 1 p @we ðCe Þ : C_ e þ S sa ma Sa n_ a P 0: e 2 @C a
ð9Þ
ð18Þ
The force-like quantity conjugate to the internal variable na is denoted Sa. Following the standard argument of Coleman and Noll [4], one obtains the elastic relation for Sp in the form
Sp ¼ 2
@we ðCe Þ : @Ce
ð19Þ
Thus the reduced dissipation inequality follows from (19) as
D :¼ a a
ð14Þ
Sa
A fundamental assumption of crystal plasticity is that plastic deformation takes place on prescribed planes. On each plane there may be more than one possible direction of plastic flow. The set of N slip systems is the set of planes and directions on which plastic deformation can occur. A slip system a 2 [1, . . . , N] is described by two orthonormal lattice vectors: the slip direction sa and the slip plane normal ma. Gurtin [7] argues that it is more meaningful in a large strain setting to use the slip rate ma, rather than the slip itself, as the primary variable for describing the plastic deformation. The single-crystal hypothesis states that the plastic distortion rate Lp is composed of the sum of the slip rates on their respective slip systems [17] p
The quasi-static balance of linear momentum in the spatial configuration is given by
Sp ¼ JFe1 rFeT :
where the elastic and plastic distortion rates are respectively defined as
Le ¼ F_ e Fe
2.3. Balance relations
ð4Þ
Consistent with the requirement that the Jacobian is positive and that plastic deformations are volume preserving, it is assumed that
L ¼ $x v
387
X
sa ma Sa n_ a P 0:
ð20Þ
a
a
Lower-case Greek superscripts denote slip-system labels. No summation convention is used for Greek superscripts. The following shorthand is used for the summation over slip systems:
4. Rate-independent and rate-dependent formulations
X
The classical theory of elastoplasticity assumes that there is a threshold value of stress (yield stress) beyond which an elastic material yields and starts deforming plastically. The elastic domain
a
¼
N X :
a¼1
ð10Þ
4.1. Rate-independent theory
388
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
E is defined as the convex region of stresses for which there results only an elastic deformation. It is assumed that there exists a set of yield functions /a(sa, Sa) (dependent on the stress and the internal variable Sa), such that /a < 0 inside the elastic domain and /a = 0 on the surface of the elastic domain. A common example of a yield function for flow on slip system a is
rate, is overdetermined. A very small rate-sensitivity parameter would require an algorithm with features developed for rate-independent systems (e.g. see Miehe and Schröder [12]).
/a ðsa ; Sa Þ ¼ jsa j Sa 6 0:
There are numerous hardening relations that govern the evolution of the flow resistance Sa. The three that will be compared in this work are outlined below.
ð21Þ
The principle of maximum plastic dissipation is the assumption that among the possible plastic configurations that satisfy the dissipation inequality (21), the actual configuration is the one that maximises the value of D [5]. The assumption of associativity between the rate measures and the normal to the yield surfaces is made, i.e.
ma
@/a ¼ ka a @s
@/a and n_ a ¼ ka ; @ðSa Þ
@/a ¼ ka sgnðsa Þ; @ sa @/a ¼ ka ; n_ a ¼ ka @ðSa Þ
ð23Þ
) n_ a ¼ jma j ¼ ka :
ð24Þ
The Karush–Kuhn–Tucker (KKT) conditions, which capture the various constraints on the yield functions and flow relations, follow as
ð25Þ
4.2. Rate-dependent theory
ma ¼ d0
a m1 js j sgnðsa Þ; Sa
ð26Þ
n_ a ¼ ma
a js j sgnðsa Þ: 1 þ 1 Sa m
ð27Þ
For a small value of the rate-sensitivity parameter m 1, the following approximations may be made for (28):
as m ! 0
1 m 1 m
þ1
1;
a js j 1; Sa
X ab h ðSÞjmb j;
a
n jm j;
ð30Þ
where S = {S1, . . . , SN}, S0 is the initial flow resistance and the hardening moduli hab describe the rate of increase of resistance on slip system a due to shear on slip system b [1]. In this theory the hardening moduli accounts for two types of hardening interactions [15]: Self-hardening, which arises from interactions between slip systems that share the same slip plane (i.e. have the same normal). Latent-hardening, which arises from interactions between slip systems that do not share the same slip plane. The coplanarity moduli vab are a useful measure of interaction between slip systems. Slip systems a and b are coplanar if mab = jma mbj = 0, and non-coplanar otherwise. The coplanarity moduli are thus
vab ¼
1 for a and b coplanar systems ðmab ¼ 0Þ for a and b non-coplanar systems ðmab – 0Þ
0
ð28Þ
ab
h ðSb Þ ¼
vab hðSb Þ þ q 1 vab hðSb Þ; |fflfflfflfflffl{zfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl}
self hardening
ð31Þ
ð29Þ
This is approximately the rate-independent result in (25). Thus, if the rate-sensitivity parameter in the constitutive Eq. (27) is very small, m 1, then the constitutive equation approximates a rateindependent model. Rate-independent models are not simple, in an algorithmic sense, to solve. This is due to the multiplicity of the slip system and that the plastic deformation rate, expressed in terms of the slip
ð32Þ
latent hardening
where q > 0 is an interaction constant and h(Sb) is a self-hardening function given by
8 < h 1 Sb a 0 b S hðS Þ ¼ : 0
for S0 6 Sb 6 S
;
ð33Þ
for Sb > S
with a, S0, S⁄ and h0 being material constants. The self-hardening function shown above was first proposed by Anand and Kothari [1]. 5.2. Fully-interactive hardening Gurtin and Reddy [8] proposed a hardening relation in the context of strain-gradient single-crystal plasticity. The special case of conventional plasticity is now considered. In the case of conventional plasticity the proposed relation has the hardening moduli dependent on the accumulated slip cbacc , rather than the slip resistance Sa. The accumulated slip on slip system b is given by
and therefore
_a
Sa jt¼0 ¼ S0 ;
b
The hardening moduli in (33) are thus given as
The rate-dependent theory arises from the regularization of the rate-independent theory. This is useful from an algorithmic standpoint, as solving a problem with multiple yield surfaces is complex [12]. A viscoplastic approximation is used to regularize the problem. A penalty term is added to the reduced dissipation inequality (21). The penalty function is chosen in the sense of a Norton–Hoff creep law [21,19]. Minimization of the resulting Lagrangian results in the following counterparts to the relations in (24):
1 m
The implicit hardening relation takes the form
ð22Þ
m a ¼ ka
/a ðsa ; Sa Þ 6 0 and ka /a ðsa ; Sa Þ ¼ 0:
5.1. Implicit hardening
S_ a ¼
where ka P 0 is the plastic multiplier. It can be shown that this assumption is equivalent to the assumption of maximum dissipation. It follows from (22) that
ka P 0;
5. Hardening relations
cbacc ¼ nb ¼
Z
t
jmb j dt:
ð34Þ
t0
As with the implicit hardening relation, the hardening moduli have contributions from the self- and latent-hardening. The coplanarity moduli (32), however, are used in the self-hardening part only. Latent hardening is intended to represent the contribution to hardening from slip systems that are not coplanar. Gurtin and Reddy argue that characterising slip-system interactions using the coplanarity moduli alone ‘‘seems too rough, because it treats all non-coplanar
389
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396 Table 1 Orientation of slip directions and normals for the double-slip crystal. Slip system
s
1
p1ffiffi ðe2 5 p1ffiffi ðe2 5
2
m p1ffiffi ð2e2 5 p1ffiffi ð2e2 5
þ 2e3 Þ 2e3 Þ
e3 Þ þ e3 Þ
where cacc :¼ ðc1acc ; . . . ; cNacc Þ. The functions SS and SL are the self- and latent-hardening contributions. They can be found by calibrating the hardening relation for the case of single-slip with the implicit relation and are given by
Fig. 1. Evolution of the flow resistance for the three calibrated hardening relations.
a
SS mslf ¼ ðS S0 Þ 1
!
1
0 1 þ S SS
SL malat ¼ qðS S0 Þ 1
h0 S
maslf
1 1þ
S S0 h0 S S
malat
ð38Þ
; ! ;
ð39Þ
with q, a, S0, S⁄ and h0 being the same material constants defined after (34). To elucidate the structure of this hardening relation an example for a crystal with two slip systems is presented in Appendix A. 5.3. Isotropic hardening In the isotropic hardening model (used by Steinmann and Stein [21], among others) the flow resistance S is governed by
SðcÞ ¼ s0 þ ðs1 s0 Þ tanh
h0 c ; s1 s0
ð40Þ
where the global accumulated slip c given by
c¼
Z
t
c_ dt ¼
t0
Fig. 2. Diagram of the indentation problem.
a
ab
i ¼ js s jm : b
ð35Þ
The hardening relation takes the form
a
a
S ðcacc Þ ¼ S0 þ SS mslf ma
þ SL mlat ; a
def
ð36Þ
ma
X
X
b
b
vab cbacc ; malat def ¼
iab cbacc ;
ð41Þ
a
s0 ¼ S0 :
ð42Þ
The relationship between the parameters in the different models will be obtained by calibration (see Section 6.1). 6. Numerical examples
where slf and lat are self- and latent-hardening measures, respectively, defined by
maslf ¼
t0
The parameters s0, s1 and h0 are material constants with s0 being the initial flow resistance. Thus
slip systems equally’’. The angle between the slip directions should also contribute to the degree with which they interact. To this end, the interaction moduli are defined by ab
Z tX jma j dt:
ð37Þ
The FEM is used to solve the system of governing equations approximately. The implementation was done using the Differential Equation Analysis Library (deal.II) ([3]). Bi-quadratic elements are used to describe the displacement variable in order to circumvent locking problems. The slip rates are solved for pointwise at the quadrature points. A predictor–corrector algorithm on the local
Fig. 3. Comparison of slip ca and flow resistance Sa for each hardening relation at the point A in the spherical indentation problem with a single slip-system crystal structure.
390
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
Fig. 4. Comparison of flow resistance Sa for each hardening relation at the point A in the spherical indentation problem with a double slip-system crystal structure.
Fig. 5. Comparison of slip ca for each hardening relation at the point A in the spherical indentation problem with a double slip-system crystal structure.
level was used alongside a global Newton–Raphson algorithm. This procedure follows that described by Steinmann and Stein [21]. The Newton–Raphson algorithm makes use of a consistent tangent in order to maintain quadratic convergence [20]. The FEM approximation that was developed was validated using the results in Steinmann and Stein [21] and de Souza et al. [5]. For all examples considered the elastic part of the free-energy was of St. Venant–Kirchhoff type and took the form
k 2 l e we Ce ¼ tr Ce I þ tr C I : 8 2
Table 2 The slip directions and normals for a face centre cubic (FCC) crystal in Miller notation. In Miller notation a normal to a plane is expressed in terms of the basis of the reciprocal lattice vectors and is denoted with round brackets, e.g. (i j k). A direction in the basis of direct lattice vectors is denoted with square brackets, e.g. [i j k]
a
sa
ma
a
sa
ma
a
sa
ma
1 2 3 4
½0 1 1
0 1
½1 0
½1 1 ½0 1 1
(1 1 1) (1 1 1) (1 1 1) 1 1Þ ð1
5 6 7 8
[1 0 1] 1 0
½1
1 1Þ ð1 1 1Þ ð1 1 1Þ ð1 1 1Þ ð1
9 10 11 12
1 0
½1 [0 1 1] 0
½1 1
1 1Þ ð1 ð1 1 1Þ ð1 1 1Þ ð1 1 1Þ
[0 1 1] ½1 0 1
[1 0 1]
ð43Þ
The rate-sensitivity parameter m from the constitutive equation for the slip rate (27) was set to m = 0.5, making the power law parabolic. The penalty parameter d0 is taken to be the initial slip rate fixed at d0 = 0.001. 6.1. Calibration of the hardening relations In order to compare the various hardening relations, the material parameters need to be calibrated. The calibration is done by considering a single slip system, where the slip resistance corre-
sponding to each hardening relation can be solved for analytically. Given a set of parameters for one hardening relation, the parameters for the other hardening relations can be found by using a nonlinear least-squares fitting algorithm. For the case of single slip, the implicit and the fully-interactive hardening relations (with a = 2 in (34)) have the same form, given by
Sðcacc Þ ¼ Ss 1 h
1 0
Ss
Ss cacc þ Ss S 0
! :
ð44Þ
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
391
Fig. 6. Two views of the distribution of the vertical component of stress (rzz) at time t = 0.2 for the spherical indentation problem with an FCC crystal structure.
Fig. 7. Two views of the distribution of slip c2 at time t = 0.2 for the spherical indentation problem with an FCC crystal structure.
Fig. 8. Comparison of flow resistance Sa for each hardening relation at point A in the spherical indentation problem with an FCC crystal structure. Note: many of the curves in (b) are near-identical.
392
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
Fig. 9. Comparison of slip ca for each hardening relation at point A in the spherical indentation problem with an FCC crystal structure.
The indenter is assumed to be perfectly bonded to the material. This results in a spherically-shaped deformation at the middle of the top of material below. Fig. 2 shows a diagram of the problem. Due to symmetry, a quarter of the crystal is modelled. A quarter of the crystal is the cube X = (0, L)3. The face z = 0 has homogeneous Dirichlet boundary conditions prescribed. Since the surface of the cube in contact with the indenter changes with time, the region on the face z = L (where the indenter acts) will change over time. The two regions on the top face are defined as
CD ðtÞ @ X :
ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 x2 þ y2 6 R2 ðL htÞ ;
z ¼ L;
ð47Þ
CN ðtÞ @ X :
ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 x2 þ y2 > R2 ðL htÞ ;
z ¼ L:
ð48Þ
Fig. 10. Diagram of the shear problem.
The parameters h0, Ss and S0 are chosen as in [1] to be
h0 ¼ 180 MPa;
Ss ¼ 148 MPa;
S0 ¼ 16 MPa:
ð45Þ
The interaction constant q = 1.4 [14,2]. For a single slip system the isotropic hardening relation takes the form
Sðcacc Þ ¼ s0 þ ðs1 s0 Þ tanh
h0 cacc : s1 s0
Here h is a parameter related to the final indentation depth. The non-homogeneous Dirichlet boundary condition on CD(t) is defined as
uðx; y; LÞ ¼
0; 0;
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R2 ðx2 þ y2 Þ þ ðL htÞ on CD ðtÞ:
ð49Þ
Homogeneous Neumann conditions on CN(t) are prescribed as
ð46Þ
rn ¼ t ¼ 0 on CN ðtÞ:
ð50Þ
The parameters h0 and s1 were found by fitting (47) to (45) with the parameters listed in (46). The value of the initial slip resistance s0 was not fitted, as it is assumed to be the same for both hardening relations, i.e. S0 = s0 = 16 MPa. The fitted parameters and a graph comparing the flow resistance evolution are shown in Fig. 1.
For each hardening relation the problem is solved using a single slip system, a double slip system and an FCC crystal with 12 slip systems. The material parameters have been calibrated to be the same for the case of a single slip system. The results from the single-slip case can thus be used to validate the calibration.
6.2. Spherical indentation
6.2.1. Single slip system Thepffiffiffi slip system normal m = e1 and slip direction s ¼ 1= 2ðe2 þ e3 Þ. As the single-slip system is used to validate the calibration of the hardening relations, it is expected that the re-
This problem concerns a rigid 0.5 L radius indenter acting on the upper surface of a single crystal over the time period t = [0, 0.2].
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
393
Fig. 11. Comparison of slip ca and flow resistance Sa for each hardening relation at the point (0.4 L, 0.4 L, 1.4 L) in the shear problem with a single slip-system crystal structure.
active formulation and the implicit formulation are the same and identical results are observed. There is a slight difference, however, with the calibrated isotropic relation. At the end of the simulation the flow resistance differs by about 2 MPa. This is most probably due to the inexactness of the fit done in the calibration. The slip, however, shows very little difference between the hardening relations. From these results it can be concluded that the calibration was successful, and therefore a comparison of the hardening relations can be made for more complicated crystal structures.
Fig. 12. Comparison of the reaction force on the top face in the direction of shear, for the three hardening relations with an FCC crystal structure.
sults for the different hardening relations should be similar. This is the case. Fig. 3 shows the slip ca and the flow resistance Sa for each of the hardening relations. For a single slip system the fully-inter-
6.2.2. Double slip system The purpose of testing the problem with a double slip system is mainly to confirm that the calibration is correct. In this case, however, identical results for the fully-interactive formulation and the implicit formulation are no longer expected. The orientation of the slip systems in the crystal are shown in Table 1. Fig. 4 shows the evolution of the flow resistance for each hardening relation. As expected, the fully-interactive formulation and the implicit formulations now give different results. The flow resistance for the first slip system seems to follow the same path for both, but the values for the second slip system differ significantly. For the implicit relation the flow resistance S2 has a value larger than S1, while for the
Fig. 13. Comparison of the ryz component of stress at time t = 0.75 s for the shear problem with an FCC crystal structure.
394
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
Fig. 14. Comparison of the slip c5 at time t = 0.75 s for the shear problem with an FCC crystal structure.
Fig. 15. Comparison of flow resistance Sa for each hardening relation at the point (0.4 L, 0.4 L, 1.4 L) in the shear problem with an FCC crystal structure. Note: many of the curves in (b) are near-identical.
fully-interactive relation it has a smaller value. This discrepancy does not carry over to the evolution of the slip, however, as shown in Fig. 5. The evolution of slip is very similar for all three hardening relations. 6.2.3. FCC crystal For the FCC crystal differences between the hardening relations are expected. The orientations of the slip systems are shown in Table 2. Figs. 6 and 7 show some views of the stress distribution and slip, solved using the fully-interactive relation, at t = 0.2. Fig. 6 shows the vertical component of the Cauchy stress r, and Fig. 7
the slip on one of the slip systems. There is a large amount of slip directly beneath the indenter, as expected. Fig. 8 shows the flow resistance for the three hardening relations. The fully-interactive formulation gives a wider flow resistance distribution. At t = 0.2 the range is from 50 MPa to 85 MPa. This is significantly different to the range of the implicit hardening relation, which is from 60 MPa to 70 MPa. Despite this large difference in flow resistance, both produce very similar results for the slip (see Fig. 9). Although there is a distinct difference between the results, they are still qualitatively similar in that they have relatively similar amounts of activity on each slip system.
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
395
Fig. 16. Comparison of slip ca for each hardening relation at the point (0.4 L, 0.4 L, 1.4 L) for the shear problem with a FCC crystal structure.
6.3. Shear problem This problem concerns a bar undergoing a shear deformation as described in [16]. The domain of the bar is X = (0, L)2 (0, 3L), L = 0.1 m, as shown in Fig. 10. Homogeneous Dirichlet boundary conditions are imposed on the z = 0 face. The displacement on the z = L face is prescribed to move in the y direction initially, then back to the original position, then in the y direction again, i.e.
8 > : ðt 23 t max Þ
if t < 13 t max if 13 tmax < t < 23 t max
ð51Þ
if t > 23 t max ;
where tmax = 2.25 is the total duration of the simulation. The remaining faces have homogeneous Neumann boundary conditions prescribed. The problem is solved for each of the three hardening relations and with two different crystal structures: the single-slip crystal and the FCC crystal. As before, the single slip system case provides validation of the calibration of the hardening relations.
6.3.1. Single slip system The slip system in this case has a normal m = e1 and slip direction s = e3. Fig. 11 shows the slip c and the flow resistance S for the three hardening relations. As before, the fully-interactive and the implicit formulations are the same and identical results are observed. A small difference in the flow resistance can be seen between the isotropic and the other two relations. Despite this difference, a good match between the slips predicted by the three hardening relations is observed. Thus it can be concluded that the calibration has been successful.
6.3.2. FCC crystal A useful measure for the shear problem is that of the reaction force on the face z = 3L. Fig. 12 shows the component of this force in the direction of shear (y direction), for each of the three hardening relations. The plastic response of the material is clear. The expansion of the yield surface due to isotropic hardening is apparent. The three hardening relations give very similar responses in this graph. Figs. 13 and 14 show the ryz component of the stress and the slip on one of the slip systems, for the fully-interactive hardening and the implicit hardening relations. The two hardening relations produce very similar results. Figs. 15 and 16 show the evolution of the flow resistance and the slip for all three hardening relations at the point (0.4L,0.4L,1.4L). The changes in the direction of shear are clear, they are the regions where the slip rate is zero in Fig. 16 and where there is an inflection in Fig. 15. The results are very similar in nature to the results from the spherical indentation problem. The flow resistance for the slip systems of the fully-interactive relation span a wider range than the implicit relation. Despite this difference, the slips are qualitatively similar for all three relations.
7. Summary and conclusion The finite element method was used to solve the equations that describe single-crystal viscoplasticity for large deformations. The algorithm used was based on that described in [21]. Three different relations for hardening were compared. The first is an isotropic relation, dependent only on the accumulated slip, as described in [21]. The second is an implicit relation, in which the rate of flow resistance is dependent on the flow resistance itself. The third hardening relation is a novel one, described in [8], in which (in
396
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
the case of conventional plasticity) the flow resistance is dependent on the slip accumulated on each slip system and takes into account fully the interaction between slip systems. In order to examine the hardening relations, two different problems were studied: a cube with the top surface being deformed by a spherical indenter and a bar undergoing a shear deformation. Of particular interest are the results for a face-centred cubic crystal structure. These show a significant difference in the evolution of the flow resistance between the fully-interactive hardening relation and the implicit relation, with the former spread over a larger range when compared to that of the implicit relation. Despite this the evolution of the slips and the reaction force for each hardening relation (including the isotropic relation) is very similar. As clearly demonstrated via the numerical examples, the fullyinteractive hardening relation can predict the same overall response as the widely-used implicit and isotropic relations. The primary advantage of the fully-interactive relation is that it can account for a wider degree of interactions between coplanar slip systems. To demonstrate that such additional complexity is required will require careful comparison with experiment. The fully-interactive hardening relation has the additional advantage, relative to the implicit one, in that that the hardening relation is dependent on the accumulated slip rather than the hardening itself. Although not emphasised in the current presentation, this form of dependence makes the numerical solution of the finite element problem more straightforward. Acknowledgement The support for this work provided by the National Research Foundation through the South African Research Chair for Computational Mechanics is gratefully acknowledged. Appendix A. Example of fully-interactive hardening relation with a double slip-system crystal structure Consider a simple two-dimensional crystal with two slip systems, the slip directions oriented at 60° to each other, i.e.
1 1 m1 ¼ pffiffiffi ð2e1 þ e2 Þ; s1 ¼ pffiffiffi ðe1 þ 2e2 Þ; 5 5 1 1 2 2 m ¼ pffiffiffi ð2e1 þ e2 Þ; s ¼ pffiffiffi ðe1 þ 2e2 Þ: 5 5
ðA:1Þ
The matrix of the coplanarity moduli takes the form
v¼
1 0 0
1
ðA:2Þ
and the matrix of the interaction moduli is
i¼
0
12 25
12 25
0
! ðA:3Þ
:
The self- and latent-hardening measures are thus
12 2 c ; 25 acc 12 1 ¼ c : 25 acc
m1slf ¼ c1acc ;
m1lat ¼
m2slf ¼ c2acc ;
m2lat
ðA:4Þ
The hardening relations for the two slip systems take the following form: S1 S2
0
!
¼ ðS S0 Þ@
0 ð1 þ S SS
ð1 q 25 12
þ
h0 S
1
0 c1acc Þ q 25 ð1 þ S SS 12
S S0 h0 12 S S 25
c
1 acc Þ
1
ð1 þ
h0 12 S 25
S S0 h0 S S
1
1
c2acc Þ þ ð1 þ q 25 Þ 12 A c
2 acc Þ
1
þ ð1 þ
q 25 Þ 12
þ
S0 S0
:
References [1] L. Anand, M. Kothari, Journal of the Mechanics and Physics of Solids 44 (1996) 525–558. [2] R.J. Asaro, A. Needleman, Acta Metallurgica 33 (6) (1985) 923–953. [3] W. Bangerth, R. Hartmann, G. Kanschat, ACM Transactions on Mathematical Software 33 (4) (2007). [4] B.D. Coleman, W. Noll, Archive for Rational Mechanics and Analysis 13 (1) (1963) 167–178. [5] E.A. de Souza Neto, D. Peric, D.R.J. Owen, Computational Methods for Plasticity: Theory and Applications, Wiley, 2011. [6] M. Gurtin, E. Fried, L. Anand, The Mechanics and Thermodynamics of Continua, Cambridge University Press, 2010. chapter 110.. [7] M.E. Gurtin, International Journal of Plasticity 24 (4) (2008) 702–725. [8] M.E. Gurtin, B.D. Reddy, Gradient single-crystal plasticity within a Mises–Hill framework based on a new formulation of self and latent hardening, 2013 (submitted for publication). [9] S.R. Kalidindi, L. Anand, C.A. Bronkhorst, Philosophical Transactions: Physical Sciences and Engineering 341 (1992) 443–447. [10] E. Kröner, Archive for Rational Mechanics and Analysis 4 (4) (1960) 273–334. [11] E.H. Lee, Journal of Applied Mechanics 36 (1969) 1–6. [12] C. Miehe, J. Schröder, International Journal for Numerical Methods in Engineering 50 (2) (2001) 273–298. [13] C. Miehe, J. Schröder, J. Schotte, Computer Methods in Applied Mechanics and Engineering 171 (3) (1999) 387–418. [14] D. Peirce, R. Asaro, A. Needleman, Acta Metallurgica 30 (6) (1982) 1087–1119. [15] D. Peirce, R.J. Asaro, A. Needleman, Acta Metallurgica 31 (12) (1983) 1951– 1976. [16] B.D. Reddy, C. Wieners, B. Wohlmuth, International Journal for Numerical Methods in Engineering 90 (6) (2012) 784–804. [17] J.R. Rice, Journal of the Mechanics and Physics of Solids 19 (6) (1971) 433–455. [18] J. Schröder, C. Miehe, Computational Materials Science 9 (1) (1997) 168–176. [19] J.C. Simo, Numerical analysis and simulation of plasticity, in: P. Ciarlet, J. Lions (Eds.), Numerical Methods for Solids (Part 3) Numerical Methods for Fluids (Part 1), Handbook of Numerical Analysis, vol. 6, Elsevier, 1998, pp. 183–499. [20] J.C. Simo, R.L. Taylor, Computer Methods in Applied Mechanics and Engineering 48 (1) (1985) 101–118. [21] P. Steinmann, E. Stein, Computer Methods in Applied Mechanics and Engineering 129 (3) (1996) 235–254.
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Finite element simulation of large-strain single-crystal viscoplasticity: An investigation of various hardening relations T.M. Povall ⇑, A.T. McBride, B.D. Reddy Centre for Research in Computational and Applied Mechanics, University of Cape Town, Private Bag, 7701 Rondebosch, South Africa
a r t i c l e
i n f o
Article history: Received 6 May 2013 Received in revised form 31 July 2013 Accepted 21 August 2013 Available online 18 September 2013 Keywords: Single crystals Large strain Plasticity Viscoplasticity Hardening Finite element method
a b s t r a c t Hardening relations describe the increase in resistance to deformation during plastic flow. Three hardening relations are compared here in the context of conventional large-strain single-crystal viscoplasticity. The first is an isotropic hardening relation. The second is a hardening relation that is expressed as an ordinary differential equation in the slip resistance. The third is a new relation, originally developed in the context of gradient crystal plasticity, in which the slip resistance is expressed explicitly in terms of the accumulated slip on each slip system. The numerical solution of the governing equations is found using the finite element method coupled with a predictor–corrector type algorithm. The features of the hardening relations are elucidated using a series of numerical benchmark problems. The parameters for the hardening relations are calibrated using a model problem. Various crystal structures are investigated, including single- and double slip, and face-centred cubic crystals. The hardening relations are compared and their relative features discussed. Ó 2013 Elsevier B.V. All rights reserved.
1. Introduction The description of plastic behaviour is characterised by (a) a yield criterion; (b) a hardening relation; and (c) a flow rule for plastic strain rates or equivalent quantities. A variety of hardening relations have been proposed for singlecrystal plasticity. Popular models include a hardening relation that is expressed as an ordinary differential equation in the flow resistance [15,9]. A new hardening relation has recently been proposed by Gurtin and Reddy [8] in the context of strain-gradient single-crystal plasticity. This relation takes into account the interactions between slip systems. The purpose of this paper is to compare the fully-interactive hardening relation [8] with an implicit relation and an isotropic relation. The fully-interactive relation was originally proposed in the context of small-strain rate-independent single-crystal gradient plasticity. The relation is extended here to the conventional largestrain viscoplastic setting. The isotropic relation is that defined by Steinmann and Stein [21] and is dependent on the sum of the accumulated slips, while the implicit relation used is that given in [9]. The finite element method (FEM) is used to solve the system of governing equations. The FEM implementation for conventional single-crystal plasticity is a well-established field with contributions from Schröder and Miehe [18], Miehe et al. [13], Steinmann and Stein [21], and Anand and Kothari [1], to mention some impor⇑ Corresponding author. Tel.: +27 (0)21 650 3817. E-mail addresses: [email protected] (T.M. Povall), [email protected] (A.T. McBride), [email protected] (B.D. Reddy). 0927-0256/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.commatsci.2013.08.043
tant examples. A predictor–corrector algorithm, similar to that described by Steinmann and Stein [21], is used in conjunction with a consistent algorithmic tangent to accurately and efficiently describe the plastic evolution. Two model problems are used to compare the hardening relations. The first involves the indentation of an elastoplastic material by a rigid spherical indenter. The second problem is the shear of a bar subjected to cyclic deformations. The various hardening relations are examined in the indentation problem for three types of crystal structures: a single slip system, a double slip system and a face-centred cubic (FCC) crystal. The shear problem examines the hardening relations using a double slip-system crystal and a FCC crystal. This paper is arranged as follows. Section 2 gives an overview of the notation, kinematics and the governing equations. Section 3 provides a short derivation of the reduced dissipation inequality. A comparison of rate-independent and rate-dependent formulations is given in Section 4. Section 5 details the three hardening relations. An explanation of the calibration process is given in Section 6.1. The results of the numerical simulation of the indentation and shear problems are discussed in Sections 6.2 and 6.3, respectively. A summary of the results and the conclusion are given in Section 7. 2. Kinematics and governing equations 2.1. Kinematics A motion of a body B is a smooth function v that assigns to each material point X at time t a point in Euclidean space, x = v(X, t). The point x is referred to as the spatial (or current) point occupied by X
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
at time t. The body B in the reference configuration at time t0 is denoted by B0. The displacement u is defined by
uðX; tÞ ¼ vðX; tÞ X:
ð1Þ
The velocity v is defined as the time rate of change of the displacement, i.e.
_ v ðX; tÞ ¼ u:
ð2Þ
The deformation gradient F is defined by
F :¼ $v ¼
@ vi ei EA ; @X A
ð3Þ
where the Jacobian J = det F is positive and ei and EA represent bases in the current and reference configurations, respectively. The Kröner–Lee decomposition multiplicatively splits the deformation gradient into an elastic and a plastic part ([10,11]), denoted Fe and Fp respectively, and is given by
F ¼ Fe Fp :
det Fe > 0 and
det Fp ¼ 1:
ð5Þ
The multiplicative split motivates the introduction of an intermediate configuration. From the elastic deformation gradient, the right elastic Cauchy–Green tensor is defined by
Ce ¼ FeT Fe :
ð6Þ
The spatial velocity gradient L is the tensor field defined by
@v i 1 ¼ ei ej ¼ Le þ Fe Lp Fe ; @xj
ð7Þ
1
1
and Lp ¼ F_ p Fp :
ð8Þ
Note that in (7) the plastic distortion rate Lp is pushed forward from the intermediate (lattice) configuration into the current configuration using Fe. 2.2. Slip systems and the single-crystal hypothesis
L ¼
X
divr þ b ¼ 0;
ð11Þ
where the Cauchy stress r is a symmetric tensor and b is the body force per unit current volume. The divergence of a second order tensor A is defined as divA = @Aij/@xjei. The resolved shear stress sa is an important measure in crystal plasticity. It is defined as the Cauchy stress resolved on slip plane a in the current configuration:
a; sa ¼ sa rm
ð12Þ
a are the slip directions and normals pushed forward where sa and m into the spatial configuration from the lattice configuration, i.e.
a ¼ Fe T ma : and m
sa ¼ Fe sa
ð13Þ
3. Free energy and dissipation inequality The dissipation inequality arises from the second law of thermodynamics [6] and is given by
X 1 J 1 Sp : C_ e þ sa ma J1 w_ P 0; 2 a
where w is the free energy and Sp is a stress measure similar to the second Piola–Kirchhoff stress, but defined in the intermediate configuration as
a
m s m :
ð15Þ
It is assumed that the free energy is dependent on both the elastic and the plastic deformation. Furthermore, it is assumed that any hardening mechanism is uncoupled from the elastic deformation. The free energy function is thus of the form
w ¼ we ðCe Þ þ wp ðna Þ;
ð16Þ
where na is a set of internal scalar plastic variables. The time derivative of the free energy is thus
@we ðCe Þ _ e X@wp ðna Þ _ a w_ ¼ n : :C þ @na @Ce a |fflfflfflffl{zfflfflfflffl}
ð17Þ
The dissipation inequality becomes, on substitution of (18),
J 1
X 1 p @we ðCe Þ : C_ e þ S sa ma Sa n_ a P 0: e 2 @C a
ð9Þ
ð18Þ
The force-like quantity conjugate to the internal variable na is denoted Sa. Following the standard argument of Coleman and Noll [4], one obtains the elastic relation for Sp in the form
Sp ¼ 2
@we ðCe Þ : @Ce
ð19Þ
Thus the reduced dissipation inequality follows from (19) as
D :¼ a a
ð14Þ
Sa
A fundamental assumption of crystal plasticity is that plastic deformation takes place on prescribed planes. On each plane there may be more than one possible direction of plastic flow. The set of N slip systems is the set of planes and directions on which plastic deformation can occur. A slip system a 2 [1, . . . , N] is described by two orthonormal lattice vectors: the slip direction sa and the slip plane normal ma. Gurtin [7] argues that it is more meaningful in a large strain setting to use the slip rate ma, rather than the slip itself, as the primary variable for describing the plastic deformation. The single-crystal hypothesis states that the plastic distortion rate Lp is composed of the sum of the slip rates on their respective slip systems [17] p
The quasi-static balance of linear momentum in the spatial configuration is given by
Sp ¼ JFe1 rFeT :
where the elastic and plastic distortion rates are respectively defined as
Le ¼ F_ e Fe
2.3. Balance relations
ð4Þ
Consistent with the requirement that the Jacobian is positive and that plastic deformations are volume preserving, it is assumed that
L ¼ $x v
387
X
sa ma Sa n_ a P 0:
ð20Þ
a
a
Lower-case Greek superscripts denote slip-system labels. No summation convention is used for Greek superscripts. The following shorthand is used for the summation over slip systems:
4. Rate-independent and rate-dependent formulations
X
The classical theory of elastoplasticity assumes that there is a threshold value of stress (yield stress) beyond which an elastic material yields and starts deforming plastically. The elastic domain
a
¼
N X :
a¼1
ð10Þ
4.1. Rate-independent theory
388
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
E is defined as the convex region of stresses for which there results only an elastic deformation. It is assumed that there exists a set of yield functions /a(sa, Sa) (dependent on the stress and the internal variable Sa), such that /a < 0 inside the elastic domain and /a = 0 on the surface of the elastic domain. A common example of a yield function for flow on slip system a is
rate, is overdetermined. A very small rate-sensitivity parameter would require an algorithm with features developed for rate-independent systems (e.g. see Miehe and Schröder [12]).
/a ðsa ; Sa Þ ¼ jsa j Sa 6 0:
There are numerous hardening relations that govern the evolution of the flow resistance Sa. The three that will be compared in this work are outlined below.
ð21Þ
The principle of maximum plastic dissipation is the assumption that among the possible plastic configurations that satisfy the dissipation inequality (21), the actual configuration is the one that maximises the value of D [5]. The assumption of associativity between the rate measures and the normal to the yield surfaces is made, i.e.
ma
@/a ¼ ka a @s
@/a and n_ a ¼ ka ; @ðSa Þ
@/a ¼ ka sgnðsa Þ; @ sa @/a ¼ ka ; n_ a ¼ ka @ðSa Þ
ð23Þ
) n_ a ¼ jma j ¼ ka :
ð24Þ
The Karush–Kuhn–Tucker (KKT) conditions, which capture the various constraints on the yield functions and flow relations, follow as
ð25Þ
4.2. Rate-dependent theory
ma ¼ d0
a m1 js j sgnðsa Þ; Sa
ð26Þ
n_ a ¼ ma
a js j sgnðsa Þ: 1 þ 1 Sa m
ð27Þ
For a small value of the rate-sensitivity parameter m 1, the following approximations may be made for (28):
as m ! 0
1 m 1 m
þ1
1;
a js j 1; Sa
X ab h ðSÞjmb j;
a
n jm j;
ð30Þ
where S = {S1, . . . , SN}, S0 is the initial flow resistance and the hardening moduli hab describe the rate of increase of resistance on slip system a due to shear on slip system b [1]. In this theory the hardening moduli accounts for two types of hardening interactions [15]: Self-hardening, which arises from interactions between slip systems that share the same slip plane (i.e. have the same normal). Latent-hardening, which arises from interactions between slip systems that do not share the same slip plane. The coplanarity moduli vab are a useful measure of interaction between slip systems. Slip systems a and b are coplanar if mab = jma mbj = 0, and non-coplanar otherwise. The coplanarity moduli are thus
vab ¼
1 for a and b coplanar systems ðmab ¼ 0Þ for a and b non-coplanar systems ðmab – 0Þ
0
ð28Þ
ab
h ðSb Þ ¼
vab hðSb Þ þ q 1 vab hðSb Þ; |fflfflfflfflffl{zfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl}
self hardening
ð31Þ
ð29Þ
This is approximately the rate-independent result in (25). Thus, if the rate-sensitivity parameter in the constitutive Eq. (27) is very small, m 1, then the constitutive equation approximates a rateindependent model. Rate-independent models are not simple, in an algorithmic sense, to solve. This is due to the multiplicity of the slip system and that the plastic deformation rate, expressed in terms of the slip
ð32Þ
latent hardening
where q > 0 is an interaction constant and h(Sb) is a self-hardening function given by
8 < h 1 Sb a 0 b S hðS Þ ¼ : 0
for S0 6 Sb 6 S
;
ð33Þ
for Sb > S
with a, S0, S⁄ and h0 being material constants. The self-hardening function shown above was first proposed by Anand and Kothari [1]. 5.2. Fully-interactive hardening Gurtin and Reddy [8] proposed a hardening relation in the context of strain-gradient single-crystal plasticity. The special case of conventional plasticity is now considered. In the case of conventional plasticity the proposed relation has the hardening moduli dependent on the accumulated slip cbacc , rather than the slip resistance Sa. The accumulated slip on slip system b is given by
and therefore
_a
Sa jt¼0 ¼ S0 ;
b
The hardening moduli in (33) are thus given as
The rate-dependent theory arises from the regularization of the rate-independent theory. This is useful from an algorithmic standpoint, as solving a problem with multiple yield surfaces is complex [12]. A viscoplastic approximation is used to regularize the problem. A penalty term is added to the reduced dissipation inequality (21). The penalty function is chosen in the sense of a Norton–Hoff creep law [21,19]. Minimization of the resulting Lagrangian results in the following counterparts to the relations in (24):
1 m
The implicit hardening relation takes the form
ð22Þ
m a ¼ ka
/a ðsa ; Sa Þ 6 0 and ka /a ðsa ; Sa Þ ¼ 0:
5.1. Implicit hardening
S_ a ¼
where ka P 0 is the plastic multiplier. It can be shown that this assumption is equivalent to the assumption of maximum dissipation. It follows from (22) that
ka P 0;
5. Hardening relations
cbacc ¼ nb ¼
Z
t
jmb j dt:
ð34Þ
t0
As with the implicit hardening relation, the hardening moduli have contributions from the self- and latent-hardening. The coplanarity moduli (32), however, are used in the self-hardening part only. Latent hardening is intended to represent the contribution to hardening from slip systems that are not coplanar. Gurtin and Reddy argue that characterising slip-system interactions using the coplanarity moduli alone ‘‘seems too rough, because it treats all non-coplanar
389
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396 Table 1 Orientation of slip directions and normals for the double-slip crystal. Slip system
s
1
p1ffiffi ðe2 5 p1ffiffi ðe2 5
2
m p1ffiffi ð2e2 5 p1ffiffi ð2e2 5
þ 2e3 Þ 2e3 Þ
e3 Þ þ e3 Þ
where cacc :¼ ðc1acc ; . . . ; cNacc Þ. The functions SS and SL are the self- and latent-hardening contributions. They can be found by calibrating the hardening relation for the case of single-slip with the implicit relation and are given by
Fig. 1. Evolution of the flow resistance for the three calibrated hardening relations.
a
SS mslf ¼ ðS S0 Þ 1
!
1
0 1 þ S SS
SL malat ¼ qðS S0 Þ 1
h0 S
maslf
1 1þ
S S0 h0 S S
malat
ð38Þ
; ! ;
ð39Þ
with q, a, S0, S⁄ and h0 being the same material constants defined after (34). To elucidate the structure of this hardening relation an example for a crystal with two slip systems is presented in Appendix A. 5.3. Isotropic hardening In the isotropic hardening model (used by Steinmann and Stein [21], among others) the flow resistance S is governed by
SðcÞ ¼ s0 þ ðs1 s0 Þ tanh
h0 c ; s1 s0
ð40Þ
where the global accumulated slip c given by
c¼
Z
t
c_ dt ¼
t0
Fig. 2. Diagram of the indentation problem.
a
ab
i ¼ js s jm : b
ð35Þ
The hardening relation takes the form
a
a
S ðcacc Þ ¼ S0 þ SS mslf ma
þ SL mlat ; a
def
ð36Þ
ma
X
X
b
b
vab cbacc ; malat def ¼
iab cbacc ;
ð41Þ
a
s0 ¼ S0 :
ð42Þ
The relationship between the parameters in the different models will be obtained by calibration (see Section 6.1). 6. Numerical examples
where slf and lat are self- and latent-hardening measures, respectively, defined by
maslf ¼
t0
The parameters s0, s1 and h0 are material constants with s0 being the initial flow resistance. Thus
slip systems equally’’. The angle between the slip directions should also contribute to the degree with which they interact. To this end, the interaction moduli are defined by ab
Z tX jma j dt:
ð37Þ
The FEM is used to solve the system of governing equations approximately. The implementation was done using the Differential Equation Analysis Library (deal.II) ([3]). Bi-quadratic elements are used to describe the displacement variable in order to circumvent locking problems. The slip rates are solved for pointwise at the quadrature points. A predictor–corrector algorithm on the local
Fig. 3. Comparison of slip ca and flow resistance Sa for each hardening relation at the point A in the spherical indentation problem with a single slip-system crystal structure.
390
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
Fig. 4. Comparison of flow resistance Sa for each hardening relation at the point A in the spherical indentation problem with a double slip-system crystal structure.
Fig. 5. Comparison of slip ca for each hardening relation at the point A in the spherical indentation problem with a double slip-system crystal structure.
level was used alongside a global Newton–Raphson algorithm. This procedure follows that described by Steinmann and Stein [21]. The Newton–Raphson algorithm makes use of a consistent tangent in order to maintain quadratic convergence [20]. The FEM approximation that was developed was validated using the results in Steinmann and Stein [21] and de Souza et al. [5]. For all examples considered the elastic part of the free-energy was of St. Venant–Kirchhoff type and took the form
k 2 l e we Ce ¼ tr Ce I þ tr C I : 8 2
Table 2 The slip directions and normals for a face centre cubic (FCC) crystal in Miller notation. In Miller notation a normal to a plane is expressed in terms of the basis of the reciprocal lattice vectors and is denoted with round brackets, e.g. (i j k). A direction in the basis of direct lattice vectors is denoted with square brackets, e.g. [i j k]
a
sa
ma
a
sa
ma
a
sa
ma
1 2 3 4
½0 1 1
0 1
½1 0
½1 1 ½0 1 1
(1 1 1) (1 1 1) (1 1 1) 1 1Þ ð1
5 6 7 8
[1 0 1] 1 0
½1
1 1Þ ð1 1 1Þ ð1 1 1Þ ð1 1 1Þ ð1
9 10 11 12
1 0
½1 [0 1 1] 0
½1 1
1 1Þ ð1 ð1 1 1Þ ð1 1 1Þ ð1 1 1Þ
[0 1 1] ½1 0 1
[1 0 1]
ð43Þ
The rate-sensitivity parameter m from the constitutive equation for the slip rate (27) was set to m = 0.5, making the power law parabolic. The penalty parameter d0 is taken to be the initial slip rate fixed at d0 = 0.001. 6.1. Calibration of the hardening relations In order to compare the various hardening relations, the material parameters need to be calibrated. The calibration is done by considering a single slip system, where the slip resistance corre-
sponding to each hardening relation can be solved for analytically. Given a set of parameters for one hardening relation, the parameters for the other hardening relations can be found by using a nonlinear least-squares fitting algorithm. For the case of single slip, the implicit and the fully-interactive hardening relations (with a = 2 in (34)) have the same form, given by
Sðcacc Þ ¼ Ss 1 h
1 0
Ss
Ss cacc þ Ss S 0
! :
ð44Þ
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
391
Fig. 6. Two views of the distribution of the vertical component of stress (rzz) at time t = 0.2 for the spherical indentation problem with an FCC crystal structure.
Fig. 7. Two views of the distribution of slip c2 at time t = 0.2 for the spherical indentation problem with an FCC crystal structure.
Fig. 8. Comparison of flow resistance Sa for each hardening relation at point A in the spherical indentation problem with an FCC crystal structure. Note: many of the curves in (b) are near-identical.
392
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
Fig. 9. Comparison of slip ca for each hardening relation at point A in the spherical indentation problem with an FCC crystal structure.
The indenter is assumed to be perfectly bonded to the material. This results in a spherically-shaped deformation at the middle of the top of material below. Fig. 2 shows a diagram of the problem. Due to symmetry, a quarter of the crystal is modelled. A quarter of the crystal is the cube X = (0, L)3. The face z = 0 has homogeneous Dirichlet boundary conditions prescribed. Since the surface of the cube in contact with the indenter changes with time, the region on the face z = L (where the indenter acts) will change over time. The two regions on the top face are defined as
CD ðtÞ @ X :
ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 x2 þ y2 6 R2 ðL htÞ ;
z ¼ L;
ð47Þ
CN ðtÞ @ X :
ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 x2 þ y2 > R2 ðL htÞ ;
z ¼ L:
ð48Þ
Fig. 10. Diagram of the shear problem.
The parameters h0, Ss and S0 are chosen as in [1] to be
h0 ¼ 180 MPa;
Ss ¼ 148 MPa;
S0 ¼ 16 MPa:
ð45Þ
The interaction constant q = 1.4 [14,2]. For a single slip system the isotropic hardening relation takes the form
Sðcacc Þ ¼ s0 þ ðs1 s0 Þ tanh
h0 cacc : s1 s0
Here h is a parameter related to the final indentation depth. The non-homogeneous Dirichlet boundary condition on CD(t) is defined as
uðx; y; LÞ ¼
0; 0;
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R2 ðx2 þ y2 Þ þ ðL htÞ on CD ðtÞ:
ð49Þ
Homogeneous Neumann conditions on CN(t) are prescribed as
ð46Þ
rn ¼ t ¼ 0 on CN ðtÞ:
ð50Þ
The parameters h0 and s1 were found by fitting (47) to (45) with the parameters listed in (46). The value of the initial slip resistance s0 was not fitted, as it is assumed to be the same for both hardening relations, i.e. S0 = s0 = 16 MPa. The fitted parameters and a graph comparing the flow resistance evolution are shown in Fig. 1.
For each hardening relation the problem is solved using a single slip system, a double slip system and an FCC crystal with 12 slip systems. The material parameters have been calibrated to be the same for the case of a single slip system. The results from the single-slip case can thus be used to validate the calibration.
6.2. Spherical indentation
6.2.1. Single slip system Thepffiffiffi slip system normal m = e1 and slip direction s ¼ 1= 2ðe2 þ e3 Þ. As the single-slip system is used to validate the calibration of the hardening relations, it is expected that the re-
This problem concerns a rigid 0.5 L radius indenter acting on the upper surface of a single crystal over the time period t = [0, 0.2].
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
393
Fig. 11. Comparison of slip ca and flow resistance Sa for each hardening relation at the point (0.4 L, 0.4 L, 1.4 L) in the shear problem with a single slip-system crystal structure.
active formulation and the implicit formulation are the same and identical results are observed. There is a slight difference, however, with the calibrated isotropic relation. At the end of the simulation the flow resistance differs by about 2 MPa. This is most probably due to the inexactness of the fit done in the calibration. The slip, however, shows very little difference between the hardening relations. From these results it can be concluded that the calibration was successful, and therefore a comparison of the hardening relations can be made for more complicated crystal structures.
Fig. 12. Comparison of the reaction force on the top face in the direction of shear, for the three hardening relations with an FCC crystal structure.
sults for the different hardening relations should be similar. This is the case. Fig. 3 shows the slip ca and the flow resistance Sa for each of the hardening relations. For a single slip system the fully-inter-
6.2.2. Double slip system The purpose of testing the problem with a double slip system is mainly to confirm that the calibration is correct. In this case, however, identical results for the fully-interactive formulation and the implicit formulation are no longer expected. The orientation of the slip systems in the crystal are shown in Table 1. Fig. 4 shows the evolution of the flow resistance for each hardening relation. As expected, the fully-interactive formulation and the implicit formulations now give different results. The flow resistance for the first slip system seems to follow the same path for both, but the values for the second slip system differ significantly. For the implicit relation the flow resistance S2 has a value larger than S1, while for the
Fig. 13. Comparison of the ryz component of stress at time t = 0.75 s for the shear problem with an FCC crystal structure.
394
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
Fig. 14. Comparison of the slip c5 at time t = 0.75 s for the shear problem with an FCC crystal structure.
Fig. 15. Comparison of flow resistance Sa for each hardening relation at the point (0.4 L, 0.4 L, 1.4 L) in the shear problem with an FCC crystal structure. Note: many of the curves in (b) are near-identical.
fully-interactive relation it has a smaller value. This discrepancy does not carry over to the evolution of the slip, however, as shown in Fig. 5. The evolution of slip is very similar for all three hardening relations. 6.2.3. FCC crystal For the FCC crystal differences between the hardening relations are expected. The orientations of the slip systems are shown in Table 2. Figs. 6 and 7 show some views of the stress distribution and slip, solved using the fully-interactive relation, at t = 0.2. Fig. 6 shows the vertical component of the Cauchy stress r, and Fig. 7
the slip on one of the slip systems. There is a large amount of slip directly beneath the indenter, as expected. Fig. 8 shows the flow resistance for the three hardening relations. The fully-interactive formulation gives a wider flow resistance distribution. At t = 0.2 the range is from 50 MPa to 85 MPa. This is significantly different to the range of the implicit hardening relation, which is from 60 MPa to 70 MPa. Despite this large difference in flow resistance, both produce very similar results for the slip (see Fig. 9). Although there is a distinct difference between the results, they are still qualitatively similar in that they have relatively similar amounts of activity on each slip system.
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
395
Fig. 16. Comparison of slip ca for each hardening relation at the point (0.4 L, 0.4 L, 1.4 L) for the shear problem with a FCC crystal structure.
6.3. Shear problem This problem concerns a bar undergoing a shear deformation as described in [16]. The domain of the bar is X = (0, L)2 (0, 3L), L = 0.1 m, as shown in Fig. 10. Homogeneous Dirichlet boundary conditions are imposed on the z = 0 face. The displacement on the z = L face is prescribed to move in the y direction initially, then back to the original position, then in the y direction again, i.e.
8 >
if t < 13 t max if 13 tmax < t < 23 t max
ð51Þ
if t > 23 t max ;
where tmax = 2.25 is the total duration of the simulation. The remaining faces have homogeneous Neumann boundary conditions prescribed. The problem is solved for each of the three hardening relations and with two different crystal structures: the single-slip crystal and the FCC crystal. As before, the single slip system case provides validation of the calibration of the hardening relations.
6.3.1. Single slip system The slip system in this case has a normal m = e1 and slip direction s = e3. Fig. 11 shows the slip c and the flow resistance S for the three hardening relations. As before, the fully-interactive and the implicit formulations are the same and identical results are observed. A small difference in the flow resistance can be seen between the isotropic and the other two relations. Despite this difference, a good match between the slips predicted by the three hardening relations is observed. Thus it can be concluded that the calibration has been successful.
6.3.2. FCC crystal A useful measure for the shear problem is that of the reaction force on the face z = 3L. Fig. 12 shows the component of this force in the direction of shear (y direction), for each of the three hardening relations. The plastic response of the material is clear. The expansion of the yield surface due to isotropic hardening is apparent. The three hardening relations give very similar responses in this graph. Figs. 13 and 14 show the ryz component of the stress and the slip on one of the slip systems, for the fully-interactive hardening and the implicit hardening relations. The two hardening relations produce very similar results. Figs. 15 and 16 show the evolution of the flow resistance and the slip for all three hardening relations at the point (0.4L,0.4L,1.4L). The changes in the direction of shear are clear, they are the regions where the slip rate is zero in Fig. 16 and where there is an inflection in Fig. 15. The results are very similar in nature to the results from the spherical indentation problem. The flow resistance for the slip systems of the fully-interactive relation span a wider range than the implicit relation. Despite this difference, the slips are qualitatively similar for all three relations.
7. Summary and conclusion The finite element method was used to solve the equations that describe single-crystal viscoplasticity for large deformations. The algorithm used was based on that described in [21]. Three different relations for hardening were compared. The first is an isotropic relation, dependent only on the accumulated slip, as described in [21]. The second is an implicit relation, in which the rate of flow resistance is dependent on the flow resistance itself. The third hardening relation is a novel one, described in [8], in which (in
396
T.M. Povall et al. / Computational Materials Science 81 (2014) 386–396
the case of conventional plasticity) the flow resistance is dependent on the slip accumulated on each slip system and takes into account fully the interaction between slip systems. In order to examine the hardening relations, two different problems were studied: a cube with the top surface being deformed by a spherical indenter and a bar undergoing a shear deformation. Of particular interest are the results for a face-centred cubic crystal structure. These show a significant difference in the evolution of the flow resistance between the fully-interactive hardening relation and the implicit relation, with the former spread over a larger range when compared to that of the implicit relation. Despite this the evolution of the slips and the reaction force for each hardening relation (including the isotropic relation) is very similar. As clearly demonstrated via the numerical examples, the fullyinteractive hardening relation can predict the same overall response as the widely-used implicit and isotropic relations. The primary advantage of the fully-interactive relation is that it can account for a wider degree of interactions between coplanar slip systems. To demonstrate that such additional complexity is required will require careful comparison with experiment. The fully-interactive hardening relation has the additional advantage, relative to the implicit one, in that that the hardening relation is dependent on the accumulated slip rather than the hardening itself. Although not emphasised in the current presentation, this form of dependence makes the numerical solution of the finite element problem more straightforward. Acknowledgement The support for this work provided by the National Research Foundation through the South African Research Chair for Computational Mechanics is gratefully acknowledged. Appendix A. Example of fully-interactive hardening relation with a double slip-system crystal structure Consider a simple two-dimensional crystal with two slip systems, the slip directions oriented at 60° to each other, i.e.
1 1 m1 ¼ pffiffiffi ð2e1 þ e2 Þ; s1 ¼ pffiffiffi ðe1 þ 2e2 Þ; 5 5 1 1 2 2 m ¼ pffiffiffi ð2e1 þ e2 Þ; s ¼ pffiffiffi ðe1 þ 2e2 Þ: 5 5
ðA:1Þ
The matrix of the coplanarity moduli takes the form
v¼
1 0 0
1
ðA:2Þ
and the matrix of the interaction moduli is
i¼
0
12 25
12 25
0
! ðA:3Þ
:
The self- and latent-hardening measures are thus
12 2 c ; 25 acc 12 1 ¼ c : 25 acc
m1slf ¼ c1acc ;
m1lat ¼
m2slf ¼ c2acc ;
m2lat
ðA:4Þ
The hardening relations for the two slip systems take the following form: S1 S2
0
!
¼ ðS S0 Þ@
0 ð1 þ S SS
ð1 q 25 12
þ
h0 S
1
0 c1acc Þ q 25 ð1 þ S SS 12
S S0 h0 12 S S 25
c
1 acc Þ
1
ð1 þ
h0 12 S 25
S S0 h0 S S
1
1
c2acc Þ þ ð1 þ q 25 Þ 12 A c
2 acc Þ
1
þ ð1 þ
q 25 Þ 12
þ
S0 S0
:
References [1] L. Anand, M. Kothari, Journal of the Mechanics and Physics of Solids 44 (1996) 525–558. [2] R.J. Asaro, A. Needleman, Acta Metallurgica 33 (6) (1985) 923–953. [3] W. Bangerth, R. Hartmann, G. Kanschat, ACM Transactions on Mathematical Software 33 (4) (2007). [4] B.D. Coleman, W. Noll, Archive for Rational Mechanics and Analysis 13 (1) (1963) 167–178. [5] E.A. de Souza Neto, D. Peric, D.R.J. Owen, Computational Methods for Plasticity: Theory and Applications, Wiley, 2011. [6] M. Gurtin, E. Fried, L. Anand, The Mechanics and Thermodynamics of Continua, Cambridge University Press, 2010. chapter 110.. [7] M.E. Gurtin, International Journal of Plasticity 24 (4) (2008) 702–725. [8] M.E. Gurtin, B.D. Reddy, Gradient single-crystal plasticity within a Mises–Hill framework based on a new formulation of self and latent hardening, 2013 (submitted for publication). [9] S.R. Kalidindi, L. Anand, C.A. Bronkhorst, Philosophical Transactions: Physical Sciences and Engineering 341 (1992) 443–447. [10] E. Kröner, Archive for Rational Mechanics and Analysis 4 (4) (1960) 273–334. [11] E.H. Lee, Journal of Applied Mechanics 36 (1969) 1–6. [12] C. Miehe, J. Schröder, International Journal for Numerical Methods in Engineering 50 (2) (2001) 273–298. [13] C. Miehe, J. Schröder, J. Schotte, Computer Methods in Applied Mechanics and Engineering 171 (3) (1999) 387–418. [14] D. Peirce, R. Asaro, A. Needleman, Acta Metallurgica 30 (6) (1982) 1087–1119. [15] D. Peirce, R.J. Asaro, A. Needleman, Acta Metallurgica 31 (12) (1983) 1951– 1976. [16] B.D. Reddy, C. Wieners, B. Wohlmuth, International Journal for Numerical Methods in Engineering 90 (6) (2012) 784–804. [17] J.R. Rice, Journal of the Mechanics and Physics of Solids 19 (6) (1971) 433–455. [18] J. Schröder, C. Miehe, Computational Materials Science 9 (1) (1997) 168–176. [19] J.C. Simo, Numerical analysis and simulation of plasticity, in: P. Ciarlet, J. Lions (Eds.), Numerical Methods for Solids (Part 3) Numerical Methods for Fluids (Part 1), Handbook of Numerical Analysis, vol. 6, Elsevier, 1998, pp. 183–499. [20] J.C. Simo, R.L. Taylor, Computer Methods in Applied Mechanics and Engineering 48 (1) (1985) 101–118. [21] P. Steinmann, E. Stein, Computer Methods in Applied Mechanics and Engineering 129 (3) (1996) 235–254.