- Email: [email protected]
A single crystal plasticity finite element formulation with embedded deformation twins
Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Contents lists available at ScienceDirect
Journal of the Mechanics and Physics of Solids journal homepage: www.elsevier.com/locate/jmps
A single crystal plasticity finite element formulation with embedded deformation twins Tao Jin a, Hashem M. Mourad a,∗, Curt A. Bronkhorst b, Irene J. Beyerlein c a
Theoretical Division, T-3, Los Alamos National Laboratory, Los Alamos, NM 87545, USA Department of Engineering Physics, University of Wisconsin, Madison, WI, 53706, USA c Department of Mechanical Engineering, Materials Department, University of California, Santa Barbara, CA, 93106, USA b
a r t i c l e
i n f o
Article history: Received 31 May 2019 Revised 23 August 2019 Accepted 10 September 2019 Available online 14 September 2019 Keywords: Deformation twinning Embedded weak discontinuity Single-crystal elasto-viscoplasticity
a b s t r a c t Deformation twinning is an important plastic deformation mechanism in some polycrystalline metals such as titanium and magnesium. In this paper, we present a novel crystal plasticity finite element framework that accounts for deformation twinning explicitly, in addition to crystallographic slip. Within this computational framework, deformation twins are treated as weak discontinuities embedded within individual finite elements, such that a jump in the velocity gradient field is introduced (via the discretized gradient operator) between the twinned and untwinned crystalline regions, taking into account compatibility and traction continuity conditions at the interface between these two regions. The deformation gradient is multiplicatively split into elastic and plastic parts in the untwinned region, as is customary in finite-deformation crystal plasticity formulations. A different multiplicative decomposition of the deformation gradient into elastic, plastic (slip), and twinning parts is adopted in the twinned region, allowing deformation twinning to be accounted for as an additional mode of plastic deformation. A stochastic model is used to predict twin nucleation at grain boundaries, and the evolution of the length and thickness of the twinned region under subsequent deformation is taken into account. The linearization of the single-crystal plasticity model, required in order to enforce traction continuity at the interface between the twinned and untwinned regions using a Newton iterative scheme, is described in detail. Simulations of the initiation and propagation of {101¯ 2} tensile twins in hexagonal close packed (HCP) titanium are presented to demonstrate the capabilities of the proposed computational framework. © 2019 Elsevier Ltd. All rights reserved.
1. Introduction The plastic deformation of metallic materials is in general mediated by the motion of dislocations, and the creation and propagation of deformation twins under mechanical loading. In order to accommodate imposed deformation, a metallic material will naturally select combinations of dislocation motion and deformation twinning. The propensity for the natural combination of the two plastic deformation mechanisms depends on the material’s crystallographic structure as well as the boundary conditions imposed upon the material to produce a given deformation field. High-symmetry materials such as copper (FCC) and tantalum (BCC) have several crystallographically close-packed atomic planes whose orientations are ∗
Corresponding author. E-mail addresses: [email protected] (H.M. Mourad), [email protected] (C.A. Bronkhorst).
https://doi.org/10.1016/j.jmps.2019.103723 0022-5096/© 2019 Elsevier Ltd. All rights reserved.
2
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
evenly spaced in three-dimensional (3D) space, allowing them to accommodate complex deformation fields via plastic slip (Clayton, 2011). Generally, the initial mobile dislocation content in such high-symmetry materials is sufficiently high, so that the resolved shear stress necessary to initiate dislocation motion is small relative to that required to nucleate deformation twins (Kubin, 2013). In addition, the Burgers vector magnitude is uniform for all systems in these crystal structures, and therefore, there is not a significant bias in resistance to dislocation motion (Bronkhorst et al., 2007; Kothari and Anand, 1998). Deformation twinning is still observed in these materials under certain loading conditions, however, generally, it is not a prominent feature. Materials such as titanium (HCP) are of low-symmetry crystallographic structure. Members of this class of materials possess a number of close-packed atomic planes, but their atomic unit cell is hexagonally cylindrical, and therefore the available slip planes/directions are not uniformly spaced in 3D. When loaded in certain directions, the unit cell does not have an adequate number of slip systems available to accommodate plastic deformation (Ardeljan et al., 2016; Beyerlein and Tomé, 2008; Gray III, 1997; Knezevic et al., 2013; Rapperport and Hartley, 1960; Salem et al., 20 03; 20 06; Song and Gray, 1995). In addition, the Burgers vector magnitude is significantly different for each of the three categories of slip systems present in HCP materials, and therefore, the resolved shear stress needed to give rise to dislocation motion also varies correspondingly (a slip system with a larger Burgers vector magnitude generally presents a greater resistance to dislocation glide; e.g. see Liu et al., 2018). For HCP materials, some estimates put these differences at several times the minimum value. The probability that deformation twinning will be required to accommodate an imposed deformation is then greater for HCP materials, and this is frequently observed experimentally (Beyerlein et al., 2010; Capolungo et al., 2009; Luan et al., 2018). Within the field of computational plasticity, deformation twinning offers unique challenges in comparison with dislocation slip. Although the physics of both mechanisms are discrete by nature, the characteristic size of deformation twins is much greater than that of dislocations (Gray, 1997; Gurao et al., 2011; Leclercq et al., 1989). Given that computational tools such as the finite element method rely upon spatial discretization of a body into individual computational elements, this disparity in characteristic size presents a challenge. Accounting for both mechanisms in full-field computational approaches such as the crystal plasticity finite element (CPFE) method therefore requires a novel approach. Although still a vigorous research area, it has become standard to represent the mechanics of dislocation motion using continuum single crystal plasticity theories (Clayton, 2011; Gurtin et al., 2010; Kubin, 2013). The characteristic length scale of dislocations still allows for an approximate homogenized representation within CPFE approaches. It still remains a challenge to adequately represent the complex interactions between dislocations and the subcell structures developed with deformation (e.g. Dequiedt et al., 2015; Hansen et al., 2010). With the characteristic length scale of deformation twins being comparable to the grain size, it is difficult to be content with a homogenized representation at the same level as for dislocations. Within the CPFE context, there is a distinct need to provide a suitable computational framework for deformation twinning, which also accounts for the evolving interfaces created. The goal of the present work is to provide this kind of differentiated representation of dislocations and deformation twins, within the context of CPFE. The process of deformation twinning is believed to involve three steps: nucleation, propagation, and growth (Beyerlein and Tomé, 2010; Cheng and Ghosh, 2017; Livescu et al., 2019; Niezgoda et al., 2014). For relatively high-purity materials, grain boundary (GB)–dislocation interaction is hypothesized to be one possible mechanism of twin nucleation (Wang et al., 20 09a; 20 09b; 2010; 2013). In this mechanism, dislocations impinging upon a GB1 interact with the latter to produce a number of partial dislocations at the impingement site (Niezgoda et al., 2014). This localized partial dislocation nucleus must be of a certain size or energy level in order for a twin to be formed within the target single crystal (Liu et al., 2018; Niezgoda et al., 2014; Wang et al., 20 09a; 20 09b; 2010; 2013). This initiation process relies upon the storage of sufficient energy within the nucleus at the GB, to allow it to propagate into the neighboring grain. In simplest form, initiation leads to the growth stage as the boundaries of the twinned region migrate in the width direction. The corresponding transformation has both shear deformation and crystallographic rotation kinematic components across the moving interface surface (Ardeljan et al., 2015; Cheng and Ghosh, 2017; Cheng et al., 2018). There have been several recent efforts to extend computational treatments of deformation twinning beyond homogenized representations. This is necessitated not only by the need to improve modeling of plasticity in general, but also by the need to integrate computational and experimental work in order to explore the science of deformation twinning more effectively. Within the domain of computational meso-mechanics for explicit representation of twinning, there has been important recent development work. Phase field techniques have recently begun to be used to examine the local evolution in structure with the inclusion of deformation twinning. The coupled processes of martensitic phase transformation together with both dislocation motion and twinning were presented by Levitas and Preston (2005) within a formal thermodynamic framework targeting NiAl, BN, C materials. The theory was formulated for large deformation; however, limited computational results were presented. A different model that represents phase transformation and deformation twinning was presented in a series of two papers (Agrawal and Dayal, 2015a; 2015b), representing a progression from one to two dimensions. This work focused on the advanced treatment of evolving interfaces associated with both transformation and twinning, represented implicitly by a continuous order parameter field. Two other publications (Clayton and Knap, 2011; 2016) firstly presented an advanced formulation for anisotropic single crystal materials along with anisotropic interface representation for twin boundaries. The
1 It is important to keep in mind that each GB is atomistically non-uniform and that dislocation impingement upon a GB surface will be discrete and non-uniform (Wang et al., 2009a, 2009b, 2010, 2013), although there is still much to be learned about the specifics.
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
3
numerical implementation was described and two-dimensional proof-of-principle calculations were also presented. In the following paper, the authors progressed to 3D polycrystal systems, while including the physics of brittle damage. Coupled twinning and damage calculations were presented to demonstrate the importance and potential of advanced meso-scale computational capability. HCP materials were the target of a recent work (Liu et al., 2018) where high-fidelity phase field representation, coupled with classical models for dislocation slip, were applied to polycrystal systems of magnesium. That work paid particular attention to nucleation at GBs, which was examined in the context of the nucleation, propagation and growth progression of interest here. HCP materials were also the focus of a series of papers (Kumar et al., 2016a; 2016b; 2019; 2015) where the spectral FFT solver was used to examine the behavior of tri-crystal configurations. Traditional finite element tools were employed by three groups to provide an explicit spatial representation of twins in 3D (Ardeljan et al., 2015; Cheng and Ghosh, 2017; Cheng et al., 2018; Chester et al., 2016; Ghosh and Cheng, 2018). This achieved significant advancements in relation to the traditional CPFE method and articulated the resolution challenges for this evolving length scale problem. The computational framework proposed here is an extension of one that has been developed previously for the representation of shear banding (Jin et al., 2019a; 2018; 2019b; Mourad et al., 2017). Recognizing that the topology of shear bands and deformation twins is similar, this computational framework is extended here to the kinematics and physics of large deformation within metallic single crystals. The current work differs from the aforementioned ones in that we propose an alternative weak-form embedding technique to explicitly represent individual deformation twins and their progression. The remainder of this paper is organized as follows. In Section 2.1, we describe the kinematics associated with the twinned and untwinned regions in the crystal and the multiplicative decomposition of the deformation gradient in each region. In Section 2.2, we discuss in detail the embedded weak discontinuity approach, which is adopted to treat deformation twinning as a heterogeneous deformation mode. In Section 2.3, we present the single-crystal hyperelastic–viscoplasticity model, and we introduce the probability-based twin nucleation model and the evolution equations of the twin growth in Section 2.4. In Section 3, we present several numerical examples to simulate the initiation and propagation of the {101¯ 2} tensile twin system in hexagonal close packed (HCP) titanium. Finally, a summary of the present paper and possible future work directions are discussed in Section 4. 2. Computational framework 2.1. Kinematics We consider a crystalline body with the reference configuration 0 undergoing the deformation ϕ over a time interval of interest T . At time t ∈ T , every particle X ∈ 0 is mapped to a corresponding point x in the deformed configuration , such that
x = ϕ ( X , t ).
(1)
The deformation gradient tensor is defined as follows:
F :=
∂ϕ = FiJ ei ⊗ E J , ∂X
(2)
where {EI } and {ei } denote the Cartesian basis vectors in the reference and current configurations, respectively. Using superposed dots to designate material time derivatives, the velocity gradient is expressed as
L :=
∂ v(x, t ) ˙ −1 = FF , ∂x
(3)
where v := u˙ is the velocity vector, defined in terms of the displacement u := x − X . Assuming that the motion of dislocations takes place on N prescribed slip systems in the crystalline lattice, the α -th α slip system (α = 1, 2, . . . , N ) is described by a pair of orthogonal unit directional vectors {t α 0 , s0 } in the reference configuα represents the slip direction. The evolution of the plastic ration, where t α represents the normal of the slip plane, and s 0 0 deformation gradient Fp is given by
F˙ p = Lp Fp ,
Lp =
α
γ˙ α Sα0 =
α
γ˙ α sα0 ⊗ t α0 ,
(4)
α where γ˙ α is the plastic slip rate on the α -th slip system. The slip plane normal t α 0 , the slip direction s0 , and the Schmid α ⊗ t α , are pushed forward according to the following relationships, tensor Sα := s 0 0 0
t α = Fe−T t α0 ,
sα = Fe sα0 ,
Sα = Fe Sα0 Fe
−1
= sα ⊗ t α ,
(5)
where tα , sα and Sα are the corresponding quantities in the deformed configuration, and Fe := FFp −1 is the elastic deformation gradient obtained via multiplicative decomposition of F, which leads to the following additive decomposition of the velocity gradient:
L = Le +
α
γ˙ α Sα , Le = F˙ e Fe −1 .
(6)
4
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Fig. 1. Multiplicative decomposition of the deformation gradient FU for the untwinned region into the elastic part FU,e and the plastic part FU,p .
Fig. 2. Multiplicative decomposition of the deformation gradient tensor FZ for the twinned region into the elastic part FZ,e , the plastic slip part FZ,p , and the twinning part FZ,tw .
After deformation twinning is initiated, plastic deformation is accommodated by both plastic slip and twinning. In this paper, we treat deformation twinning as a heterogeneous deformation mode, and represent it using the embedded weakdiscontinuity approach detailed in Section 2.2. We use (•)U and (•)Z to represent quantities associated with the untwinned and twinned crystal regions, respectively. Thus, in the untwinned region where plastic deformation is accommodated by dislocation glide on the prescribed slip systems, the deformation gradient FU is decomposed, as illustrated in Fig. 1, in accordance with the classical multiplicative decomposition (Kröner and Teodosiu, 1974; Lee, 1969; Lee and Liu, 1967):
FU = FU,e FU,p ,
(7)
where FU,p characterizes the local plastic deformation due to dislocation motion and FU,e represents the elastic stretching and rotation of the crystal lattice. The plastic part of the deformation gradient, FU,p , evolves according to Eq. (4); i.e.
F˙ U,p = LU,p FU,p ,
LU,p =
α
γ˙ α Sα0 =
α
γ˙ α sα0 ⊗ t α0 .
(8)
In the twinned region, deformation twinning is included as an additional mode of plastic deformation, as illustrated in Fig. 2. As a result, the multiplicative decomposition of the deformation gradient takes the following form (Addessio et al., 2017; Feng et al., 2018; Tjahjanto et al., 2008):
FZ = FZ,e FZ,p FZ,tw .
(9) FZ,tw ,
The twinning part of the deformation gradient, which maps the reference configuration onto the twinning configuration (the first relaxed intermediate configuration shown in Fig. 2), is given by η
FZ,tw = 1 + γtw mη ⊗ nη ,
(10) η
where 1 is the second-order identity tensor, and γtw , mη and nη are, respectively, the magnitude of the twinning shear
strain, the shearing direction, and the unit normal to the habit plane, characteristic of the (η-th) twin system under consideration. Furthermore, in Eq. (9), the elastic stretching and rotation of the crystalline lattice inside the twin is represented by FZ,e . Plastic deformation due to dislocation motion inside the twinned region is characterized by FZ,p , which evolves according to
F˙ Z,p = LZ,p FZ,p ,
LZ,p =
α
α
γ˙ α Sˆ 0 =
α
α
γ˙ α sˆ α0 ⊗ tˆ 0 .
(11) α
α
Here, LZ,p represents the plastic velocity gradient in the twinned region. The slip plane normal tˆ 0 , the slip direction sˆ 0 and α α ˆ 0 := sˆ α ˆ the Schmid tensor S 0 ⊗ t 0 , corresponding to the α -th slip system in the twinned region of the crystal, are related to
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
5
their counterparts in the untwinned region via α
α
tˆ 0 = R t α0 ,
sˆ 0 = R sα0 ,
α
ˆ 0 = R Sα0 RT , S
(12)
where R is an orthogonal matrix representing the rotation between the lattice orientations in the untwinned and twinned regions, given by (Kalidindi, 1998; Van Houtte, 1978):
R = 2 n 0 ⊗ n 0 − 1,
(13)
in terms of the unit vector n0 , normal to the twinning plane in the reference configuration. It is important to note that the tensorial quantities used to describe the deformation in the untwinned and twinned regions are not completely independent as they need to satisfy compatibility conditions at the interface formed between these two regions. In addition, the stresses in the crystal must satisfy traction continuity (equilibrium) conditions at the parent– twin interface. In Section 2.2, we discuss these constraints in more detail and describe the methods used to enforce them within the present computational framework, where deformation twins are represented as embedded weak discontinuities. 2.2. Embedded weak discontinuity finite element formulation Representation of deformation twinning in a computational setting presents unique challenges compared to dislocation slip, mainly owing to the large disparity in length scales associated with these two plastic deformation mechanisms. While the length scale associated with dislocations allows for an approximate homogenized representation in CPFE simulations, it is difficult to be content with a volume fraction-based, pseudo-slip representation of twinning, in view of the fact that its characteristic length scale is comparable to the grain size. One limitation of such homogenized treatments of twinning, for instance, is related to their inability to resolve the local stress state at the twin–parent grain interface, driving the growth of an existing twin. The need for a computational approach which can address such limitations has led to the development of CPFE approaches where twins are resolved explicitly using conventional finite element methods coupled with single crystal plasticity models (Ardeljan et al., 2015). Such an approach, however, requires remeshing of every grain containing a twin, immediately after the twin is formed and at every subsequent time step while the twin is growing, which can render CPFE simulations prohibitively expensive, especially since a highly refined mesh is often needed to discretize twinned regions due to their typical high aspect ratio. To address the aforementioned computational challenges, we adopt an embedded weak-discontinuity finite element formulation, based on a generalization of the Hu–Washizu variational principle to dynamic problems, where the velocity, velocity gradient, and stress are treated as three independent fields. This allows twinning to be treated as a heterogeneous deformation mode, represented within an individual finite element as a discontinuity in the velocity gradient field (i.e. a weak discontinuity). In the present section, we show how the assumed velocity gradient, which appears in the Hu–Washizu variational statement of the problem, is designed to admit such weak discontinuities, such that twins can be embedded within the finite element grid (mesh) instead of being resolved explicitly. This leads to a sub-grid finite element formulation, capable of representing features of the solution (in this case, weak discontinuities representing deformation twinning) that are smaller than the size of the finite element grid, which obviates the need for highly refined meshes and expensive remeshing operations. Given that twinning tends to be more pronounced in dynamic problems involving high-strain rate loading conditions, and since the twin propagation velocity is comparable to the shear wave speed in the material, requiring extremely small time steps in order to resolve the relevant time scales, the formulation includes inertial effects and adopts an explicit time stepping scheme. 2.2.1. Variational form Introducing the assumed velocity gradient L and the assumed stress S, and treating these as independent field variables along with the velocity v, the three-field variational form is expressed as follows,
δ =
δv · (ρ v˙ − ∇ · S ) d +
+
δ S : ( L − L ) d +
h
δ L : ( σ − S ) d
δv · (S · n¯ − h¯ ) d = 0,
(14)
revealing the corresponding Euler–Lagrange equations,
ρ v˙ − ∇ · S = 0 S = σ ( L, q ) L=L S · n¯ = h¯
(balance of linear momentum), (constitutive relation), (compatibility condition), (Neumann boundary condition).
(15)
Here, h¯ is the applied traction acting upon the subset h of the domain’s boundary, which has an outward unit normal vector denoted by n¯ . Recall that L = ∇ v is the velocity gradient, as defined in Eq. (3) and note that σ = σ (L, q ) is the stress obtained through the constitutive relations, where q denotes internal state variables. It is also noted that Eq. (14) is the same updated Lagrangian weak form used for the treatment of shear banding in previous works (Jin et al., 2019a; 2018; 2019b; Mourad et al., 2017). While it is possible to develop an alternative total Lagrangian weak form, this additional development
6
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
effort is not undertaken herein since the two Lagrangian formulations are essentially identical in terms of accuracy and computational performance (see Belytschko et al., 2013, pp. 147–148). 2.2.2. Interpolation scheme The deformed configuration is partitioned into ne disjoint element domains, e , with boundaries e , and the following interpolation scheme for the fields v, L, L, and S, is adopted:
vi = NiA d˙ A ,
(16)
Li j = B¯ i jA d˙ A ,
(17)
˜ i jA d˙ A , Li j = B
(18)
Si j = Si jA sA ,
(19)
where d˙ A denotes a velocity degree of freedom,2 NiA is the corresponding nodal shape function, sA is a stress degree of ¯ operator in freedom, Si jA is a stress interpolation function, and summation is implied over repeated indices. Using the B the discretization of L is equivalent to selective reduced integration of volumetric and shear strain terms (see Mourad et al., 2017), and thus alleviates volumetric locking arising in the presence of well-developed (isochoric) plastic flow, and mitigates parasitic shear effects manifested under bending-dominated conditions. It is noted that the assumed velocity gradient L is ˜ , described in detail below. discretized using a different operator, B 2.2.3. Finite element equations Substituting Eqs. (16)–(19) into (14) and writing every integral over the global domain (or the subset h of its boundary) as a sum of integrals over each element e (or subset he of the element’s boundary), we obtain
δ d˙ A
MAB d¨ B + fAint − fAext + δ sA HAB d˙ B = 0,
(20)
for all arbitrary variations δ d˙ A and δ sA , i.e., MAB d¨ B + fAint = fAext ,
(21)
HAB d˙ B = 0,
(22)
where MAB =
ne
e
e=1
HAB =
ne
fAint =
e
ne e=1
fAext =
e=1
e
ne
he
e=1
ρ NiA NiB d,
(23)
¯ i jB − B ˜ i jB d, Si jA B
(24)
˜ i jA σi j d + HBA sB , B
(25)
NiA h¯ i d .
(26)
As shown by Simo and Hughes (1986), ensuring that Eq. (22) holds for any admissible velocity field requires (L − L ) to be orthogonal to the admissible stress field S, i.e. H = 0. This orthogonality condition is expressed as follows:
e
¯ −B ˜ d = 0, B
(27)
assuming a piecewise constant stress field. As a result, Eq. (25) is simplified, and the internal force is given by fAint =
ne e=1
e
˜ i jA σi j d. B
(28)
It should be noted here that σ ij is the stress obtained through the constitutive equations as a function of the assumed velocity gradient L. 2 Here we adopt the notation of Belytschko et al. (1988); Fish and Belytschko (1988). This notation reflects our computer implementation where degrees of freedom are stored in one-dimensional arrays. For example, individual components of nodal velocity vectors are stored in the one-dimensional array d˙ , and the index of this array is denoted by a capital letter, A or B, and encapsulates both the node number and the component number (i.e. direction in space) corresponding to each nodal degree of freedom.
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
7
Fig. 3. In the present approach, a deformation twin is modeled as a weak discontinuity embedded within individual finite elements; i.e. a strain jump is introduced across the interface between the embedded twin RZ and the surrounding parent grain (untwinned region) RU .
2.2.4. Assumed velocity gradient field In the embedded weak discontinuity approach, the velocity (and displacement) fields remain continuous, while discontinuities are introduced into the velocity gradient (strain rate) field to describe strain jumps across the interface between the embedded twin and the surrounding material (untwinned region), as shown schematically in Fig. 3. The present embedded discontinuity approach is classified as a sub-grid method as it allows a twinned region of width b ≤ h to be embedded within a finite element with characteristic size h. The jump in the velocity gradient field assumes the following mathematical form (Fish and Belytschko, 1988):
[[L]] = g m ⊗ n = g Ti j ei ⊗ e j ,
(29)
where g is the magnitude of the jump, n is the unit vector normal to the twin plane, and m is a unit vector such that T = m ⊗ n represents the discontinuous deformation mode. Noting that L: T represents the projection of the velocity gradient onto the discontinuous mode, the assumed velocity gradient is written in component form as follows:
Li j = Li j + gTi j =
δik δ jl + α (x )Ti j Tkl Lkl ,
(30)
αZ
where δ ij is the Kronecker delta, and α (x) is a scalar parameter such that α (x ) = in the embedded twin and α (x ) = −α U ˜ is given by: in the untwinned region. Accordingly, in view of Eqs. (17)–(18), the discontinuous operator B
δik δ jl + α Z Ti j Tkl B¯ klA in the embedded twin, ˜ i jA = B ˜U B = δik δ jl − α U Ti j Tkl B¯ klA in the untwinned region, i jA
˜ Zi jA = B
(31)
¯ and unit normal to the twinning plane n, and the unknown parameters m, in terms of its known continuous counterpart B α Z and α U , which are determined as described below. Integrals over the element domain e , such as the one in Eq. (28), are approximated via Gaussian quadrature as follows:
e
f ( x ) d ≈
Nq
wq J ( ξ q ) f ( ξ q ),
(32)
q=1
where Nq is the number of integration points, ξ q and wq are the coordinates of the q-th integration point in e and its corresponding weight coefficient, respectively, and J is the Jacobian determinant of the isoparametric mapping. After the twin is embedded within an element, the integrand f(ξ q ) assumes two different values at each integration point, fZ (ξ q ) and fU (ξ q ), representing the embedded twin and the surrounding untwinned material, respectively. The following modified quadrature rule, introduced by Belytschko et al. (1988) and Fish and Belytschko (1988) , is adopted in this case to evaluate the total integral:
e
f ( x ) d ≈
Nq
wq J ( ξ q )
q=1
b Z f (ξ q ) + h
1−
b h
f U (ξ q ) ,
(33)
√ where b is the width of the embedded twin, and h = A is the characteristic size of the element e , whose area is denoted by A. Based on the above modified quadrature rule, the orthogonality condition (27) can be expressed as:
e
¯ −B ˜ d ≈ B
Nq q=1
wq J (ξ q ) B¯ (ξ q ) −
b Z ˜ (ξ q ) + B h
1−
b ˜ U (ξ q ) B h
= 0.
(34)
Substituting (31) into (34), we obtain the following expression for the parameter α Z in terms of α U :
α = Z
h −1 b
αU.
(35)
8
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Similarly, based on Eq. (28), the element internal force vector is written as follows:
fAint
e
=
˜ i jA σi j d ≈ B
e
Nq
wq J ( ξ q )
q=1
b ˜ Zi jA (ξ q ) σiZj (ξ q ) + B h
1−
b U ˜U B i jA (ξ q ) σi j (ξ q ) , h
(36)
where σ Z and σ U denote the Cauchy stress in the twin and the untwinned region, respectively. Noting that the same slipbased mechanisms govern the evolution of plastic deformation in these two regions, the stresses σ Z and σ U are computed using the same single-crystal plasticity model described in Section 2.3. It is important to note, however, that the slip systems α {Sˆ 0 } and {Sα0 }, on which dislocation motion takes place in the twinned and untwinned regions, respectively, have different orientations, and the plastic deformation rates, LZ,p and LU,p , in these two regions are different as well; cf. Eqs. (8) and (11). Equilibrium at the interface between the embedded twin and the surrounding untwinned material requires σ Z and σ U to satisfy the following traction-continuity condition:
n ·
Z Z,e Z,p Z,p σ F , F , L − σ U FU,e , FU,p , LU,p d = 0
(37)
where d is an element of area on the twin boundary in the deformed configuration. The above nonlinear equation is solved for m and α U using the Newton iterative scheme detailed in Appendix A. 2.2.5. Time integration The time interval of interest, T , is split into increments [t n , t n+1 ] of size t := t n+1 − t n . Quantities evaluated at the beginning and end of a typical time increment are denoted, respectively, by (•)n and (• )n+1 . Each time step, the nodal accelerations d¨ n+1 are obtained by solving the dynamic finite element equations (21) explicitly: −1
ˆ d¨ n+1 = M
int f ext n+1 − f n+1 ,
(38)
ˆ for improved computational efficiency, and then the replacing the consistent mass matrix M by its lumped counterpart M central difference scheme is used to update the nodal solution variables as follows:
1 d¨ n t 2 , 2
(39)
1 d¨ n + d¨ n+1 t. 2
(40)
dn+1 = dn + d˙ n t + d˙ n+1 = d˙ n +
The global internal force vector f int n+1 in Eq. (38) is computed via standard finite element assembly procedures (Hughes, 20 0 0) from its elemental counterpart given by Eq. (36) in terms of the updated stresses σ Zn+1 inside the twin and σ Un+1 in the untwinned material. These stresses are in turn computed using the single-crystal elasto-viscoplasticity model detailed in several published works (Anand, 2004; Kalidindi et al., 1992; Kalidindi, 1998) and outlined in Section 2.3. This constitutive model requires, as input, a deformation gradient that consists of elastic and plastic (slip) parts only. The appropriate deformation gradient tensors are obtained as follows: Recasting Eq. (3) as F˙ = LF, and using the mid-point exponential map integrator, we obtain
Fn+1 = exp[ L]Fn ≈ [1 + L]Fn ,
(41)
where
L := tLn+ 12 = t
∂ vn+ 12 ∂ u ∂ (un+1 − un ) = = . ∂ xn+ 12 ∂ xn+ 12 ∂ xn+ 12
(42)
As noted above, in the present finite element formulation, the stress must be evaluated using the constitutive model as a function of the assumed velocity gradient L. Therefore, based on (41), we write
Fn+1 ≈ [1 + L]Fn
⇒
Fn+1 =
FZn+1 ≈ 1 + LZ FZn FU n+1
≈ 1 + L
U
FU n
in the embedded twin, in the untwinned region,
(43)
where, in light of Eq. (30),
Z Li j = δik δ jl + α Z Ti j Tkl Lkl in the embedded twin, L i j = LUi j = δik δ jl − α U Ti j Tkl Lkl in the untwinned region.
(44)
The stress in the untwinned material, σ U , is computed using the single-crystal plasticity model as a function of FU , n+1 n+1 which can be decomposed into elastic and plastic (slip) parts according to Eq. (7). On the other hand, the stress within the twin, σ Zn+1 , is computed using the same material model as a function of
Fˆ Zn+1 := FZn+1 FZ,tw
−1
,
which can be decomposed multiplicatively into elastic and plastic (slip) parts according to Eq. (9).
(45)
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
9
2.3. Constitutive modeling Based on the kinematic relationships described in Section 2.1 and the assumption that thermally activated slip is the dominant mechanism for plastic deformation in the twinned and untwinned regions alike, the stress in these two regions is computed using the same single-crystal elasto-viscoplasticity model. This constitutive model is outlined below, and a more detailed description thereof can be found in the literature (Anand, 2004; Kalidindi et al., 1992; Kalidindi, 1998). Although the model includes temperature as a state variable, all simulations presented herein are conducted assuming isothermal conditions. For convenience, we employ notation that does not differentiate between quantities associated with the twinned or untwinned regions in the following brief description of the material model (Section 2.3.1) and its algorithmic implementation (Section 2.3.2). These differences are highlighted in Section 2.3.3, and the material Jacobian is discussed in Section 2.3.4.
2.3.1. Single-crystal elasto-viscoplastic model The theoretical foundations of single-crystal elasto-viscoplasticity can be traced back to the work of Asaro and Needleman (1985) Asaro and Rice (1977), Asaro (1983), Hill and Rice (1972), and Rice (1971). Let Cc represent the elasticity tensor in the crystal basis that linearly depends on the temperature θ as follows:
Cci jkl (θ ) = C0i jkl + Mi jkl θ ,
(46)
where C0i jkl is the elasticity tensor at θ = 0 K, and Mi jkl are constant coefficients. Let Q represent the orthogonal tensor that rotates the crystal basis to the global basis. Then, the fourth-order elasticity tensor in the global basis is given by
Ci jkl = Qim Q jn Qkp Qlq Ccmnpq .
(47)
The free energy is expressed as follows:
1 2
ψ ( Ee ) = Ee : C ( θ ) : Ee ,
(48)
where Ee is the elastic Green–Lagrange strain tensor defined in the relaxed (second) intermediate configuration as:
Ee =
1 e 1 ( C − 1 ) = Fe T Fe − 1 . 2 2
(49)
The second Piola–Kirchhoff (2nd P-K) stress Se in the (second) intermediate configuration, which is work-conjugate to the elastic strain tensor Ee , is calculated as:
Se = C : Ee .
(50)
Recall that the plastic deformation rate Lp in Eq. (4) is based on the Schmid tensor and the plastic slip rate. The flow rule for the slip rate γ˙ α on the α -th slip system is expressed as (Bronkhorst et al., 2007; Busso, 1990):
F γ˙ α = γ˙ 0 exp − 0 1 − kθ
|τ α | − sαρ μ/μ0 sαl μ/μ0
p q
sign(τ α ),
(51)
where k is the Boltzmann constant, γ˙ 0 is a constant reference rate of shearing, sα ρ is the deformation resistance due to evolving dislocation density, sα is the constant intrinsic lattice resistance, the Macaulay brackets • designate the ramp l function, i.e. x := max(x, 0 ), ∀x ∈ R, the shear modulus ratio μ(θ )/μ0 encapsulates the temperature sensitivity, and τ α is the resolved stress on the α -th slip system, defined as:
τ α = (Ce Se ) : Sα0 .
(52)
The evolution of the deformation resistance sα ρ on slip system α is described as:
s˙ αρ =
hαβ |γ˙ β |,
hαβ = r + (1 − r )δαβ hβ ,
β
β
β
ss − sρ hβ = h0 β β ss − s0
(53)
where r is the latent-hardening parameter (Peirce et al., 1982), hβ is the self-hardening rate including both the dislocation β generation and annihilation effects (Acharya and Beaudoin, 20 0 0), and ss is the saturation stress parameter that depends on the temperature and shear rate as follows (Kocks, 1976): β
β
ss = ss0
|γ˙ β | γ˙ 0
kAθ .
(54)
10
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
2.3.2. Algorithmic implementation We use the compact notation t = tn and τ = tn+1 to represent, respectively, the beginning and the end of a typical time interval [tn , tn+1 ]. Letting γ α = γ˙ α (τ ) t represent the increment of the plastic shearing on the α -th slip system, and using Eq. (4) we write:
F (τ ) = exp t p
α
γ˙ α (τ )Sα0 Fp (t ) = exp
Therefore,
Fp (τ )−1 = Fp (t )−1 exp −
α
α
γ α Sα0 Fp (t ) ≈
1+
α
γ α Sα0 Fp (t ).
γ α Sα0 ≈ Fp (t )−1 1 − γ α Sα0 .
(55)
(56)
α
Based on the multiplicative decomposition, F(τ ) = Fe (τ )Fp (τ ), we have:
F ( τ ) = F ( τ )F ( τ ) e
p
−1
≈ F(τ )F (t ) p
−1
1−
α
γ α Sα0 .
(57)
Assume that the deformation increment F = F(τ )F(t )−1 in the time step is purely elastic, we define the trial elastic deformation at the end of the time step as follows:
Fe trial (τ ) = FFe (t ) = FF(t )Fp (t )−1 = F(τ )Fp (t )−1 .
(58)
As a result, the relationship between Fe (τ ) and Fe trial (τ ) can be expressed as:
Fe (τ ) ≈ Fe trial (τ ) 1 −
α
γ α Sα0 .
(59)
Based on Eq. (59), the right Cauchy–Green tensor is expressed as:
C (τ ) = F (τ ) F (τ ) ≈ e
e
T e
1−
=
1−
α
α
γ α (Sα0 )T Fe trial (τ )T Fe trial (τ ) 1 −
α
γ α Sα 0
γ α (Sα0 )T Ce trial (τ ) 1 − γ α Sα0
≈ Ce trial (τ ) −
α
α
γ α Ce trial (τ )Sα0 + Sα0 T Ce trial (τ ) ,
(60)
where the second-order term O (( γ α )2 ) is neglected. Therefore, the elastic strain is written as:
Ee (τ ) ≈ Ee trial (τ ) −
α
γ α sym Ce trial (τ )Sα0 ,
(61)
and the 2nd P-K stress Se (τ ) at the end of the time step is given by:
Se (τ ) = C : Ee (τ ) ≈ C : Ee trial (τ ) − = Se trial (τ ) −
α
α
γ α sym Ce trial (τ )Sα0
γ α C α ( τ ) ,
(62)
where Cα (τ ) = C : sym[Ce trial (τ )Sα ]. 0 α In a single crystal, the prescribed slip systems in the reference configuration {t α 0 , s0 }, {α = 1, 2, . . . , N }, are timeindependent. In the time interval [t, τ ], the deformation gradients at the beginning and the end of the time step, F(t) and F(τ ), are known. In a finite element simulation, the role of the constitutive update subroutine is to compute the variables {Fp (τ ), sαρ (τ ), Se (τ )} at the end of the time step, starting from their known values {Fp (t ), sαρ (t ), Se (t )} at the beginning of the time step. This update procedure is governed by a set of implicit equations, summarized as follows, solved using a two-level iterative procedure detailed by Anand (2004) and Kalidindi et al. (1992):
Se (τ ) = Se trial (τ ) − α γ α Cα (τ ), sαρ (τ ) = sαρ (t ) + β hαβ (τ )| γ β |, γ α = γ˙ α (τ ) t.
(63)
2.3.3. Differences between twinned and untwinned regions For conciseness, a notation that does not differentiate between quantities associated with the twinned and untwinned regions has been employed in the foregoing overview of the constitutive equations and stress update procedure. For the sake of clarity, several important points are emphasized here:
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
11
1. The tensors Ee and Se are defined on the intermediate configuration in the untwinned region, and on the second intermediate configuration in the twinned region; cf. Figs. 1 and 2. α 2. When the resolved shear stress τ α is computed according to Eq. (52), the Schmid tensors {Sα } and {Sˆ 0 } are used in 0 the untwinned and twinned regions, respectively. 3. The elasticity tensors CU and CZ corresponding to the untwinned and twinned regions, respectively, are related as follows (in the global basis):
CZi jkl = Rim R jn Rkp Rlq CU mnpq ,
(64)
where R is the orthogonal matrix, defined in Eq. (13), which relates the lattice orientations in the twinned and untwinned regions of the crystal. 4. In the stress update procedure, FU = FU,e FU,p is provided as input to compute the stress σ U in the untwinned region, while Fˆ Z := FZ (FZ,tw )−1 = FZ,e FZ,p is provided as input to compute the stress σ Z inside the twin. 2.3.4. Material Jacobian In order to enforce the traction continuity condition (37), the Newton scheme described in Appendix A requires the calculation of the material Jacobian W (τ ) based on the constitutive equations summarized in the present section,
W (τ ) =
∂ σ (τ ) , ∂ F (τ )
(65)
where σ is the Cauchy stress tensor obtained by pushing forward the 2nd P-K stress Se as follows:
σ=
1 T Fe Se Fe . detFe
(66)
The derivation of the material Jacobian W is presented in Appendix B. Note that the Jacobian derived by Balasubramanian (1998) and Dai (1997), is defined as follows:
¯ := W
∂σ ∂σ ≈ , ∂E ∂U
based on the approximation
E = lnU ≈ U − 1, where U is the right stretch tensor obtained via polar decomposition of the deformation gradient, F = RU. In Appendix B, we derive the expression for the Jacobian defined in Eq. (65), which is required in the embedded-discontinuity finite element formulation. 2.4. Twin nucleation, propagation and growth models The twin initiation criteria adopted in this paper is based on a probabilistic twin nucleation model detailed in (Beyerlein and Tomé, 2010; Niezgoda et al., 2014). One of the core assumptions of this probability-based model is that twin nucleation occurs when some number of grain boundary dislocations undergo stress-driven transformations that then coalesce into a single stable nucleus. The cumulative total number of such transformations N(τ ) on the stress interval [0, τ ] is treated as a random variable and described by a Poisson process. The probability of observing N = k transformations on the stress interval [0, τ ] is expressed as:
P (N (τ ) = k ) =
(τ )k k!
e−(τ ) .
(67)
In Eq. (67), (τ ) represents the cumulative rate of the grain boundary dislocation transformations, defined as:
(τ1 , ) =
τ1 + τ1
λ ( τ ) dτ ,
(68)
where is the stress interval, and λ(τ ) represents the probability of the grain boundary dislocation transformations. In order to generate a stable twin nucleus, it is assumed that at least k∗ transformation events are required and this reorganization is a stress driven process. Let the random variable S represent the critical nucleation strength, or the stress required to transform an appropriate number of grain boundary dislocations, then,
P (τ > S ) = P (N (τ ) ≥ k ) = 1 − ∗
∗ k −1
P ( N ( τ ) = k ).
(69)
k=0
Based on the assumption that k∗ = 1 and all dislocation events generate a stable twin nucleus, we have
P (τ > S ) = P (N (τ ) ≥ 1 ) = 1 − e−(τ ) .
(70)
12
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
In order to decide the specific form of the cumulative rate of the grain boundary dislocation transformations (τ ), the corresponding probability λ(τ ) is assumed to follow the Weibull distribution, so that λ(τ ) → 0 as τ → 0, preventing the case of non-zero transformation rate at zero stress. As a result,
(τ ) =
τ β τc
,
(71)
where τ c is the Weibull scale parameter and can be interpreted as a characteristic stress for nucleation, and β is the Weibull shape parameter. Both τ c and β are treated as material constants and do not depend on temperature or loading conditions. In the context of single-crystal elasto-viscoplasticity presented in Section 2.3, τ is the resolved shear stress associated with α ˆ 0 under consideration. In the finite element simulation of twin propagation using the embedded weak the twin system S discontinuities introduced in Section 2.2, once the twin nucleation probability associated with an element, calculated according to Eq. (70), exceeds the threshold value P0 , that is,
P (τ > S ) = 1 − e−(τ ) > P0 ,
(72)
we assume that the twin has nucleated inside this element. In the finite element simulation, all the state variables, including the twin nucleation probability, are evaluated at the Gauss quadrature points. In the current implementation, the twin nucleation probability within an element is taken as the highest value calculated at the Gauss points inside the element. An alternative approach is to evaluate Eq. (72) at the element centroid. After twin nucleation is detected at a grain boundary, the generated twin starts to grow, and eventually propagates into neighboring elements under subsequent deformation. Let ltw represent the embedded twin length inside individual elements, and recall that b represents the embedded twin thickness. We adopt the following expressions for the twin propagation rate l˙tw and the twin thickness growth rate b˙ tw :
−Gtw l˙tw = vl exp 1− kθ
τ p q sl
,
−Gtw b˙ tw = vb exp 1− kθ
τ p q sl
,
(73)
where vl is the reference propagation speed, and vb is the reference thickness growth speed. During a typical time interval [tn , tn+1 ], the embedded twin length and the twin thickness at the end of the time step are calculated as: tw ˙tw lntw +1 = ln + ln+1 t,
bn+1 = bn + b˙ tw n+1 t.
(74)
Once the twin length ltw inside an element is larger than the length of the line segment le bounded by the intersection points between the twin plane and the element boundary, that is,
l tw ≥ l e ,
(75)
the twin is allowed to propagate into the neighboring element. We summarize the computational framework, including the twin initiation and propagation scheme, in Algorithm 1 , where we use the following terminology and notation. The grain boundary set, denoted by Egb , is the set of elements located on the grain boundaries between different neighboring crystals. The embedded twin set, denoted by Etw , is the set of e , is the set of neighboring elements elements containing the embedded twin. The neighboring element set, denoted by Enb e e that share a common edge with an individual element . The sets of Egb and Enb are known at the start of the finite element simulation from the initial arrangement of multi-crystals and the spatial discretization, and the set of Etw is empty before twin nucleation is detected. 3. Numerical results In this section, we demonstrate the capability of the proposed computational framework through several numerical examples in two dimensions (2D). Specifically, we focus on simulating the initiation and propagation of the {101¯ 2} tensile twin system in Titanium. In the first numerical example, we construct a simple tri-crystal structure, and assume that a fully developed {101¯ 2} tensile twin is embedded in the center of the middle crystal (Grain C) from the beginning of the simulation. We aim to verify the embedded weak discontinuity formulation by comparing the corresponding numerical results with the ones obtained from the conventional finite element simulation that explicitly resolves the embedded tensile twin. In the second numerical example, we use the same tri-crystal structure. However, instead of assuming that the location of the twin is known a priori, we rely on the probability-based twin nucleation criterion presented in Section 2.4 to detect the initiation of the twin. Moreover, instead of assuming that the twin is fully developed from the beginning, we allow the twin to evolve (propagate and grow) according to the growth rates expressed in Eq. (73). In the third numerical example, we simulate the initiation and propagation of the {101¯ 2} tensile twin in a more general poly-crystalline structure. The material parameters used for HCP titanium are given in Table 1. A constant temperature of 295 K is assumed everywhere in the domain in all the simulations. The simulations for all the aforementioned problems are conducted in two-dimensions. In reality, however, structures of crystal lattices, and their slip and twin systems, are three-dimensional. Therefore, spatial rotations are necessary to simplify the {101¯ 2} twin propagation that is intrinsically three-dimensional, to a two-dimensional setup. Let Tc represent a secondorder tensor defined in the orthogonal basis {eci } associated with the crystal lattice, and T represent the counterpart defined
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
13
Algorithm 1: Main steps performed within a single time increment [tn , tn+1 ].
1
2 3
• Loop over elements in the discretized domain h : if e ∈ Egb (current element is located at the grain boundary) then Calculate the twin nucleation probability Pi according to Eq. (70) at i-th Gauss point in e ; if max {Pi } > P0 (P0 is the threshold value) then Add e to the embedded twin set Etw ; Compute the twin plane normal in the current configuration via push-forward of its (known) counterpart in the undeformed configuration: n =
4 5
6
7 8 9 10 11
12
F −T n n+1 0
;
F −T n n+1 0
Store the coordinate of the element centroid as x0 ; Obtain the twin propagation trajectory {x|(x − x0 ) · n = 0}; if e ∈ Etw (current element has embedded twin) then , α Z and α U by solving the traction continuity condition (37) using the Newton scheme detailed in Calculate m Appendix A; Z and σ U , in the twinned and untwinned regions, respectively; Update stresses σ Compute the contribution of e to the internal force fAint e , using Eq. (36); Calculate the twin propagation rate l˙itw and the thickness growth rate b˙ tw according to Eq. (73) at i-th Gauss point; i Update the twin length litw and the twin thickness btw inside the element according to Eq. (74) at i-th Gauss point; i Calculate the length of the line segment l e bounded by the intersection points between the twin propagation trajectory and the boundary of the element, {x|(x − x0 ) · n = 0} ∩ ∂ e ; if max {litw } ≥ l e then
• Loop over the neighboring elements enb ∈ Enb of the current element e :
if enb ∈ / Etw and x|x ∈ enb ∩ {x|(x − x0 ) · n = 0} = ∅ (the neighbor element enb intersects with the twin
13
14 15
propagation trajectory) then Add enb to the embedded twin set Etw ; else U in the untwinned material at all Gauss points in e ; Update stress σ Compute the contribution of e to the internal force fAint e , using Eq. (36) with b = 0;
Table 1 Material parameters used for HCP titanium. Parameter
Value
Unit
Parameter
Value
Unit
ρ θ
4507 298.0 177.59 191.34 87.57 45.01 52.23 −52.0 −36.7 13.3 −32.8 −19.7 1.0 × 107 250.0
kg/m3 K GPa GPa GPa GPa GPa GPa/K GPa/K GPa/K GPa/K GPa/K s−1 MPa
F0 p q r h0 s0 A ss0
3.0 × 10−19 0.6 1.4 1.4 225.0 30.0 1.0 × 10−19 400.0 11.1 200.0 0.85 3190.0 10.0 2.76 × 10−19
J — — — MPa MPa J MPa — MPa — m/s m/s J
C011 , C022 C033 C012 , C013 , C023 C044 C055 , C066 M11 , M22 , M33 M12 , M13 , M23 , M44 M55 , M66 γ˙ 0 sl
β τc
P0 vl vb Gtw
in the global Cartesian coordinate basis {ei }. The two tensors T and Tc are related through a rotation matrix Q in the following form:
T = Q Tc QT . The rotation matrix Q rotates the crystal basis c
e = Qe .
(76)
{eci }
to the global basis {ei } according to the following relationship:
(77)
The rotation represented by the matrix Q can be decomposed into a (non-unique) sequence of three simple rotations using Euler angles {φ , θ , ω}. The three simple counter-clockwise rotations adopted in this paper are enumerated as follows (Kalidindi et al., 1992):
14
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Fig. 4. The tri-crystal problem, where the fully developed {101¯ 2} tensile twin is embedded in the center of Grain C at the beginning of the simulation: (a) the boundary conditions, (b) the velocity profile applied at the top surface (corresponding to a shear strain rate of 7.5 × 104 s−1 ).
Fig. 5. Rotation of the crystal lattice so that, initially (at t = 0), the plane spanned by the normal and shear directions (n and m) of the {101¯ 2} tensile twin system coincides with the e1 –e2 plane (i.e. the x–y plane) in the global coordinate system.
(1) Given an orthogonal basis {x, y, z} aligned with {ec1 , ec2 , ec3 } associated with the crystal lattice, rotate {x, y, z} by the first Euler angle φ about the z-axis to get {x(1) , y(1) , z}. This is represented by:
Q1 =
cos φ −sin φ 0
sin φ cos φ 0
0 0 . 1
(2) Rotate {x(1) , y(1) , z} by the second Euler angle θ about the x(1) -axis to get {x(1) , y(2) , z(1) } so that z(1) is aligned with the out-of-plane direction e3 . This is represented by:
Q2 =
1 0 0
0 cos θ −sin θ
0 sin θ . cos θ
(3) Rotate {x(1) , y(2) , z(1) } by the third Euler angle ω about the z(1) -axis to get {x(2) , y(3) , z(1) } so that it is aligned with {e1 , e2 , e3 } of the global basis. This is represented by:
Q3 =
cos ω −sin ω 0
sin ω cos ω 0
0 0 . 1
Combining the three simple rotations described above, the rotation matrix Q is expressed as:
Q = Q3 Q2 Q1 =
cos φ cos ω − sin φ cos θ sin ω −cos φ sin ω − sin φ cos θ cos ω sin φ sin θ
sin φ cos ω + cos φ cos θ sin ω −sin φ sin ω + cos φ cos θ cos ω −cos φ sin θ
sin θ sin ω sin θ cos ω . cos θ
(78)
3.1. Tri-crystal with a pre-existing twin We construct the tri-crystal problem shown in Fig. 4, where three grains with different orientations undergo a simple shear in the x-direction. The tri-crystal has a total width of 200 μm and a height of 50 μm. Grains A and B have a width of 50 μm each. Fig. 4(a) shows the boundary conditions used in the simulation, where the horizontal velocity profile shown in Fig. 4(b) is applied to the top surface of the tri-crystal sample, corresponding to a shear strain rate of 7.5 × 104 s−1 . Based
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
15
Fig. 6. With the conventional FEM, (a) the {101¯ 2} tensile twin is explicitly resolved using a refined finite element mesh. With the embedded weak discontinuity approach, (b) a coarser mesh is used in the parent grains/untwinned region, and (c) the twin is represented by weak discontinuities embedded within a subset of the elements (highlighted in blue). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 7. Contours of the shear stress σ 12 (unit: MPa) at very early stages in the simulation, showing as the stress wave (a) begins propagating from the top surface, (b) reaches the twinned region, (c) reaches the bottom surface, and then (d) bounces back and starts propagating upward.
on the rotation matrix defined in Eq. (78), the Euler angles {φ = 210◦ , θ = 90◦ , ω = 42.472◦ } are used to rotate the crystal lattice associated with Grain C (the middle grain), as shown in Fig. 5, so that the plane spanned by the normal and shear directions {n, m} of the {101¯ 2} tensile twin system in Grain C coincides with the x–y plane in the global coordinate system. Thus, before the tri-crystal deforms (at t = 0), the normal to the twin plane points in the y-direction, i.e. n = [0, 1, 0]T , and the shear direction points in the x-direction, i.e. m = [1, 0, 0]T , in the global coordinate system. Grains A and B are both assigned the Euler angles {φ = 210◦ , θ = 90◦ , ω = 27.472◦ }, so that two 15◦ tilt grain boundaries are formed (see Fig. 4). With the aid of these crystal rotations, the problem is cast in a two-dimensional (2D) setting. We assume that the {101¯ 2} tensile twin has fully developed inside Grain C from the beginning, as shown in Fig. 4, and the location of the twin is known a priori (in the middle of Grain C). Based on the transformation matrix between the lattice orientations in the untwinned and the twinned regions, as shown in Eq. (13), the Euler angles of the crystal lattice inside the embedded twin is {φ = 30◦ , θ = 90◦ , ω = −42.472◦ }. In order to verify the proposed computational framework, we use two different approaches in the numerical simulations. In the first approach, we use a refined mesh to explicitly resolve the tensile twin, as shown in Fig. 6(a). The numerical results obtained using this approach are treated as the reference solution. In the second approach, we adopt a coarser mesh and use the proposed computational framework to treat the tensile twin as embedded weak discontinuities inside Grain C, as shown in Figs. 6(b) and (c). Two cases are considered. In the first case, the thickness b of the twinned region is 4.5 μm, and thus the ratio between the twin thickness and the element characteristic length is close to unity (b/h ≈ 1.0). Fig. 7 shows the propagation of the
16
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Fig. 8. Distributions of the third Euler angles ω (unit: degree) in Grains A and B, and in the twinned and untwinned regions of Grain C. Results obtained using the embedded weak discontinuity appraoch with b/h ≈ 1.0 (b and c) are in very good agreement with those obtained using the conventional FEM (a).
Fig. 9. Comparison of the conventional FEM and the embedded weak discontinuity simulations: (a) evolution of the third Euler angle ω in the twinned and untwinned regions, and (b) the load–displacement curves.
shear stress σ 12 in very early stages of the simulation conducted using the proposed computational framework. As expected, once the velocity profile shown in Fig. 4(b) is applied to the top surface, the stress wave propagates downward until it reaches the bottom surface, and then bounces back. Fig. 8 shows the distributions of the third Euler angle ω in the different grains at the end (t = 10−6 s) of the simulations conducted using the conventional finite element method and the embedded weak discontinuity approach. The lattice rotations in the twinned and untwinned regions obtained using the two numerical approaches are in very good agreement, which is confirmed by the evolution history of the third Euler angle ω reported in Fig. 9(a). Moreover, the load–displacement plots generated using the embedded weak discontinuity and conventional FEM results, are indistinguishable as can be seen from Fig. 9(b). Fig. 10 further examines other state variables, including the cumulative slip and the stress components inside the twinned and untwinned regions, obtained using the two numerical approaches. Again, the numerical results obtained using the embedded weak discontinuity approach, shown in Fig. 10(b), (c), (e), (f), (h), (i), (k) and (l), are in good agreement with the corresponding results obtained using the conventional FEM with a refined mesh that explicitly resolves the twin, shown in Fig. 10(a), (d), (g) and (j). Note that among the three stress components σ 11 , σ 22 and σ 12 , only the first one (σ 11 ) is different in the twinned region (Fig 10(f)) compared to the untwinned region (Fig. 10(e)). The other two stress components (σ 22 and σ 12 ) have the same magnitudes in the twinned and untwinned regions; see Fig. 10(h), (i), (k) and (l). This is due to the enforcement of the traction continuity condition (37), along the (horizontal) interface between the twinned and untwinned regions. In the second case, the thickness b of the fully developed twin is 1.5 μm, and the ratio between the twin thickness and the element characteristic length is b/h ≈ 1/3. Figs. 11 and 12 compare the results obtained using the conventional FEM and
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
17
Fig. 10. Distributions of the cumulative slip on all slip systems (a–c) and the in-plane stress components (d–l, unit: MPa) obtained using the conventional FEM and the embedded weak discontinuity approach with b/h ≈ 1.0. The results obtained using the proposed framework (b, c, e, f, h, i, k, l) are in very good agreement with the corresponding results obtained using the conventional FEM (a, d, g, j).
Fig. 11. Distributions of the 3rd Euler angles ω (unit: degree) in Grains A, B, and the twinned and untwinned regions of Grain C. The numerical results obtained using the proposed framework with b/h ≈ 1/3 (c–f) are in good agreement with the corresponding results obtained using conventional FEM (a–b). Note that the band of elements highlighted in blue in (e–f) is about three times wider than the twinned region embedded within these elements. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
18
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Fig. 12. Distributions of the cumulative slip on all slip systems (a–c) and the in-plane stress components (d–l, unit: MPa) obtained using the conventional FEM and the embedded weak discontinuity approach with b/h ≈ 1/3. The numerical results obtained using the proposed framework (b, c, e, f, h, i, k, l) are in good agreement with the corresponding results obtained using conventional FEM (a, d, g, j). Note that the band of elements highlighted in blue in (e–f) is about three times wider than the twinned region embedded within these elements. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 13. Twin nucleation probability distribution in the elements adjacent to the grain boundaries (a) at t = 0.60 × 10−6 s, and (b) at t = 0.646 × 10−6 s.
the embedded weak discontinuity approach. Note that in these two figures, distributions of state variables in the twinned region are plotted over the whole area of the elements in which the twin is embedded, even though the twinned region only occupies one third of the area of those elements. As in the first case, we observe good agreement between the results of the embedded weak discontinuity approach and the conventional FEM simulation for the naive tri-crystal problem, in which the fully developed twin is embedded into the mesh at the beginning of the simulation. The above numerical simulations are used as a simple test to verify the implementation of the proposed computational framework.
3.2. Twin initiation and propagation in tri-crystals In this numerical example, we use the same geometry and boundary conditions of the previous tri-crystal problem, shown in Fig. 4. However, in this case, we adopt the probability-based twin nucleation model, introduced in Section 2.4, to detect the location of twin initiation at the grain boundaries. Moreover, instead of assuming that the twin is fully developed from the beginning, we allow the twin to evolve according to the evolution Eq. (73). That is, after twin initiation is detected at one grain boundary, the twin thickness and length both grow under subsequent deformation (b and ltw are functions of time t), and the twin propagates into neighboring elements once inequality (75) is satisfied. Fig. 13 shows the distribution of the twin nucleation probability at the two grain boundaries. Once initiation of the {101¯ 2} tensile twin is detected at one of the grain boundaries (at t = 0.323 × 10−6 s), the twin starts to grow in thickness and propagate along the twinning plane into the neighboring grain. Fig. 14 shows the twin propagation and the growth of the twin thickness b(t) during the deformation process. After the twin is initiated at the left grain boundary, it starts to
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
19
Fig. 14. Growth and propagation of the {101¯ 2} tensile twin during the deformation process in the tri-crystal problem. Note that the characteristic element size in this problem is h ≈ 5.7 μm.
Fig. 15. Stress components σ 11 and σ 12 (unit: MPa) in the untwinned (a, c) and twinned (b, d) regions at t = 0.660 × 10−6 s.
thicken while at the same time propagating to the right. This propagation process continues until the twin reaches the right grain boundary. Fig. 15 shows the distributions of the stress components σ 11 and σ 12 inside the untwinned and twinned regions at t = 0.660 × 10−6 s, approximately when the twin reaches the mid-point between the left and right grain boundaries. During this process, the traction continuity condition n · [[σ ]] = 0 is enforced along the horizontal interface between the twin and
20
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Fig. 16. Evolution of the stress components σ 11 and σ 12 (unit: MPa) inside the untwinned (left) and twinned (right) regions after the twin is fully developed: The σ 11 component exhibits stress relaxation in the embedded twin (b, f, j) but remains almost constant in the untwinned region (a, e, i); meanwhile, the σ 12 component remains equilibrated ([[σ12 ]] = 0) between the twinnned and untwinned material, due to proper enforcement of the traction continuity condition n · [[σ ]] = 0.
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
21
Fig. 17. Poly-crystal problem: (a) the geometry and boundary conditions, (b) the velocity profile applied at the top surface, corresponding to a strain rate of 15.0 × 105 s−1 , and (c,d) Finite element mesh and value of the third Euler angle ω in each grain at t = 0.
the surrounding material. This can be confirmed, for example, by comparing Fig. 15(c) and (d), and observing that σ 12 has the same value in the twinned region and the surrounding untwinned material (i.e. the jump [[σ12 ]] = 0). Conversely, since the interface is horizontal (i.e. n points in the x2 -direction), traction continuity does not impose any constraint on the σ 11 component, and thus the jump [[σ11 ]] = 0, as can be seen by comparing Fig. 15(a) and (b). Fig. 16 shows the distributions of the same two stress components inside the untwinned and twinned regions at three different points in time, and shows that the same observations hold true after the conclusion of the twin propagation stage. It is also observed by comparing Fig. 15(a) to (b), and Fig. 16(a) to (b), that initially, σ 11 is larger inside the twin than in the surrounding untwinned material. However, stress relaxation takes place very shortly thereafter inside the twinned region, as can be seen by comparing Fig. 16(b), (f), and (j). This stress relaxation is due to dislocation slip-based plastic deformation, which is concentrated inside the twinned region as can be seen by comparing Fig. 10(b) to (c), and Fig. 12(b) to (c). The short delay associated with this effect can be attributed to the rate-dependent nature of the plastic response, represented by the single crystal constitutive model described in Section 2.3.1.
22
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Fig. 18. Twin nucleation probability distributions in the elements adjacent to the grain boundaries at different time steps. Two twin nucleation events are detected. Results obtained using a mesh size h = 2 μm (a–c) are compared to those obtained using a coarser mesh with h = 4 μm (d–f).
3.3. Twin initiation and propagation in poly-crystal In this example, the numerical simulation is performed on a 100 × 100 μm2 poly-crystalline structure that undergoes uniaxial tension, as shown in Fig 17. Figs. 17(a) and (b), respectively, show the boundary conditions and the velocity profile applied at the top surface, corresponding to a strain rate of 15 × 105 s−1 . A constant temperature of 295 K is assumed everywhere in the domain. Two different meshes, with the average mesh size h = 4 μm and h = 2 μm, are used to discretize the poly-crystalline domain, as shown in Fig. 17(c) and (d). Each grain is assigned a different third Euler angle ω in order to form tilt grain boundaries. Fig. 18 compares the distributions of twin nucleation probability at the grain boundaries obtained from simulations using different mesh sizes. Fig. 18(a), (b) and (c) show the results obtained at different time steps from the simulation with the mesh size h = 2 μm, and Figs. 18(d), (e) and (f) show the results obtained at different time steps from the simulation with the mesh size h = 4 μm. Fig. 18 shows that two twin nucleation events are detected at the grain boundaries. Due to the influence of the mesh size h on the magnitude of the stress calculated at the grain boundary, the second twin nucleation event is detected at slightly different time values, depending on the mesh size used in the simulation. However, by comparing Fig. 18(b) with (e) and Fig. 18(c) with (f), we can see that the first and second twin nucleation sites are mesh-independent. After twin initiation is detected, the two sets of twins start to propagate and grow in thickness. Fig. 19 shows the propagation paths of the tensile twins initiated at the grain boundaries. In this figure, the dark gray area represents elements adjacent to the grain boundaries, and the red area represents the elements containing the embedded twin. It is worth emphasizing that the size of the elements highlighted in red does not represent the thickness of the embedded twin. Once initiated from the grain boundaries, the two tensile twins start to propagate into the middle grain, and the propagation continues until they reach another grain boundary. We also compare the twin propagation paths obtained from the simulations using different mesh sizes. Fig. 19(a), (b) and (c) are obtained from the simulation using the mesh size h = 2 μm, and
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
23
Fig. 19. Propagation paths of the two twins initiated at the grain boundaries. The red areas consist of elements containing an embedded twin, and the dark gray areas comprise elements that are adjacent to the grain boundaries. Results obtained using a mesh size h = 2 μm (a–c) are compared to those obtained using the coarser mesh with h = 4 μm (d–f). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 19(d), (e) and (f) are obtained from the simulation using the mesh size h = 4 μm. It can be seen from these figures that the twins propagate along the same paths, regardless of the mesh size used in each simulation. During the twin propagation process, the thickness of the deformation twin evolves according to Eq. (73), which accounts for the temperature and the resolved shear stress on the twin system under investigation. Fig. 20 shows contour plots of the twin thickness at various points in time during the deformation process, where the field plotted in color is the thickness b of the embedded weak discontinuity contained within each element. It is noted that the patch of elements containing the embedded twin has serrated (zig-zag) boundaries due to the finite element spatial discretization. Also note that this patch of elements has a different thickness depending on the mesh size used. However, it can be seen from Fig. 20 that the average twin thickness is approximately 0.3 μm, irrespective of the (significantly larger) mesh sizes used in the two cases considered (2 μm and 4 μm). In other words, by comparing the results obtained from simulations using two different mesh sizes, we conclude that the evolution of the embedded twin thickness does not depend on the underlying mesh size. It is also emphasized that, consistent with experimental observations, the embedded twin follows a straight propagation path which intersects (i.e. passes through) the set of finite elements forming the aforementioned serrated patch. Fig. 21 compares the distributions of the σ 22 stress in the twinned and untwinned crystal regions after the deformation twin is fully developed. Numerical results obtained from simulations using different mesh sizes are compared, as done previously. It is clear from Fig. 21, that the middle grain has a relatively smooth stress distribution inside and outside the embedded twin, and that there are no visible stress concentrations around the patch of elements containing the embedded twin. 4. Summary In this paper, we present a new computational approach to incorporate deformation twinning as a plastic deformation mechanism, in addition to crystallographic slip, into an explicit finite element formulation. Different deformation gradients are adopted for the untwinned and twinned crystal regions based on the underlying kinematics, which is a departure from the conventional isodeformation gradient assumption that is widely used in the literature. More specifically, in the
24
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Fig. 20. Evolution of the twin thickness b within the elements that contain the embedded twin. Numerical results obtained using the mesh size h = 2 μm (a-c) are compared to those obtained using the mesh size h = 4 μm (d–f). Note that it is the value of b within each element, not the size h of the element itself, that determines the thickness of the embedded twin.
untwinned crystal region, plastic deformation is dominated by plastic slip and the deformation gradient is multiplicatively decomposed into elastic and plastic parts. This is consistent with large deformation crystal plasticity formulations. In the twinned crystal region, deformation twinning is considered as an additional plastic deformation mechanism, and the associated deformation gradient is multiplicatively decomposed into elastic, plastic and twinning parts. The two deformation gradients are related through the compatibility condition and the traction continuity condition at the interface between the untwinned and twinned regions. In the finite element formulation, the deformation twin is modeled as a heterogeneous deformation mode, and treated as weak discontinuities embedded within individual elements. Multiple twins are allowed to nucleate at the grain boundaries within a poly-crystalline structure by using a probability-based twin nucleation model. Once initiated at the grain boundaries, twin propagation and growth are simulated according to the corresponding evolution equations under subsequent deformation. Several numerical examples are provided to verify the implementation, and to demonstrate the capabilities of the proposed crystal plasticity finite element framework with regard to the accurate, mesh-independent representation of the deformation twinning process. The computational framework presented in this paper is only the first step toward incorporating deformation twinning in general three-dimensional (3D) finite element simulations of crystal plasticity. Future improvements can be achieved at least through the following avenues. The single-crystal elasto-viscoplasticity model adopted in this paper is based on a purely mechanical theory and assumes isothermal conditions, i.e. although the temperature is treated as a state variable, it is not allowed to evolve during the deformation process. In order to depict the underlying physics related with deformation twinning more accurately, a constitutive model that adequately represents thermo-mechanical effects will be adopted in future work. The current constitutive model can easily be replaced by such a coupled thermo-mechanical model without significant changes to the structure of the proposed computational framework itself. For all the simulations presented in this work, deformation twinning is simplified as a two-dimensional deformation process via spatial rotations of the crystal lattices under investigation. However, this plastic deformation mechanism, like crystallographic slip, is intrinsically threedimensional. Therefore, more sophisticated numerical methods, especially with regard to the evolution of the twin topology, are required to ensure the propagation continuity of the deformation twin in a general three-dimensional setup. Addition-
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
25
Fig. 21. Distributions of the stress component σ 22 in (a,b) the untwinned crystal regions and (c,d) the twinned crystal regions after the twin is fully developed in the middle grain (t = 60.0 × 10−9 s) using different mesh sizes.
ally, the propagation and growth of deformation twins are highly complex multi-physical processes that involve different length-scales. Combined theoretical and experimental research efforts are required to arrive at an accurate description of this physical phenomenon. Acknowledgments Funding for this work was provided by the Laboratory Directed Research and Development (LDRD) program of Los Alamos National Laboratory under project number 20150431ER, by the DoD/DOE Joint Munitions Program (JMP), and by the Advanced Simulation and Computing–Physics and Engineering Models (ASC–PEM) program. The authors wish to acknowledge fruitful discussions about the present work with Hansohl Cho, Biao Feng, Marko Knezevic, Veronica Livescu, and Jason R. Mayeur. Appendix A. Newton scheme for the traction continuity equation The traction continuity condition (37) can be recast as:
n · [[σ ]] = 0,
(A.1)
where [[σ [[ represents the jump of the stress field at the interface between the embedded twin and the surrounding material. Recall that during the time interval [tn , tn+1 ], t = tn+1 − tn represents the time increment, (•)n represents the variable
26
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
at time tn that is already known, and (• )n+1 represents the counterpart at time tn+1 that needs to be updated. The above nonlinear equation is to be solved using a Newton scheme. Let (•)(k) and (• )(k+1 ) represent the quantities associated with the k-th and (k + 1 )-th iterations. Based on Eq. (43), we have
δ FUn+1 = δ LU FUn ,
δ FZn+1 = δ LZ FZn .
(A.2)
By using the mid-point rule and the relationship between α Z and α U shown in Eq. (35), we recast the rate form between L and L, as shown in Eq. (30), into the incremental form:
LZi j = δik δ jl + (h/b − 1 )α U Ti j Tkl Lkl in the embedded twin, L i j = LUi j = δik δ jl − α U Ti j Tkl Lkl in the untwinned region.
(A.3)
ˆ , defined as follows, Introducing the variables αˆ and m
αˆ = α U Tkl Lkl , mˆ = αˆ m,
(A.4)
we rewrite Eq. (A.3) in component form:
LUi j = Li j − mˆ i n j ,
LZi j = Li j +
or in tensorial form:
LU = L − mˆ ⊗ n,
L Z = L +
h ˆ in j , −1 m b
(A.5)
h ˆ ⊗ n. −1 m b
Substituting into Eq. (A.2), we obtain
(A.6)
δ FUn+1 = δ LU FUn = −δ mˆ ⊗ n FUn , δ FZn+1 = δ LZ FZn =
h −1 b
δ mˆ ⊗ n FZn .
(A.7)
The Cauchy stress is linearized, yielding the following expression for the stress at the end of the (k + 1 )-th iteration at t = tn+1 : +1 ) ) ) ) σ n(k+1 ≈ σ n(k+1 + δσ n(k+1 = σ n(k+1 +
) ∂ σ n(k+1 ) ) : δ Fn+1 = σ n(k+1 + W n(k+1 : δ Fn+1 , ∂ Fn+1
(A.8)
where W represents the material Jacobian associated with the single-crystal elasto-viscoplasticity material model of Section 2.3. The derivation of the Jacobian W is presented in Appendix B. Using the linearization of the stress shown in Eq. (A.8), the traction continuity condition (A.1) is linearized as follows: +1 ) n · [[σ ]]n(k+1 =n ·
σZ − σU
(k+1) n+1
≈n ·
σ Z + δσ Z
(k )
n+1
−
σ U + δσ U
(k ) n+1
= 0.
(A.9)
For notational convenience, we drop the subscript indicating the time step and rearrange the above equation into the following form,
n·
δσ Z − δσ U
(k )
= −n ·
σZ − σU
(k )
(A.10)
According to Eqs. (A.7) and (A.8), we have
δσ Z = W Z : δ FZ = W Z : and
h −1 b
δ mˆ ⊗ n FZn ,
(A.11)
δσ U = W U : δ FU = −W U : δ mˆ ⊗ n FUn .
(A.12)
Substituting Eqs. (A.11) and (A.12) into (A.10), we obtain (in component form):
ni
h Z U − 1 WiaZ jk Fn,pk + WiaU jk Fn,pk b
( k )
ˆ j = −ni n pδm
σiaZ − σiaU
(k )
.
(A.13)
This can be recast into the following compact form:
ˆ j = g(ak ) . za(kj ) δ m
(A.14) (k )
where the (k + 1 )-th iteration residual ga
g(ak ) = −ni za(kj ) = ni
U (k )
σiaZ − σia
(k )
and tangent matrix za j are given by
,
h Z U − 1 WiaZ jk Fn,pk + WiaU jk Fn,pk b
(A.15)
( k ) n p.
(A.16)
ˆ . Finally, we recover the values of α U and The linearized system of Eq. (A.14) is used to solve iteratively for the vector m ˆ with the aid of Eq. (A.4). m from the converged solution for m
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
27
Appendix B. Derivation of the material Jacobian Our goal here is to derive the Jacobian of the single-crystal elasto-viscoplasticity model, defined as follows:
W (τ ) =
∂ σ (τ ) . ∂ F (τ )
This material Jacobian is used to solve the nonlinear traction continuity condition, using the Newton scheme described in Appendix A. Quantities at the beginning of the time step (•)t are known. Dropping the time indicator τ associated with quantities at the end of the time step for convenience, and expressing all tensors in indicial notation, we have:
∂σi j ∂ 1 e e e = F S F ∂ Fkl ∂ Fkl detFe im mn jn e 1 e ∂ Fjne ∂ Fim e e 1 e ∂ Smn e e e e e e ∂ = S F + F F + F S + F S F , im im mn jn detFe ∂ Fkl mn jn ∂ Fkl jn im mn ∂ Fkl ∂ Fkl detFe
Wi jkl =
and
1 ∂ 1 ∂ Fe e −T = − F : . ∂ Fkl detFe detFe ∂ Fkl
(B.1)
(B.2)
Letting
T =
∂ Fe ( τ ) , ∂ F (τ )
Q=
∂ Se ( τ ) , ∂ F (τ )
(B.3)
we have
Wi jkl =
1 e e e e e e e e e e −1 Timkl Smn Fjn + Fim Qmnkl Fjn + Fim Smn T jnkl − Fim Smn Fjn Fqp T pqkl . e detF
(B.4)
• Calculation of T :
Ti jkl =
∂ Fiej ∂ p = Fim Fmn (t )−1 δn j − γ α Sα0n j ∂ Fkl ∂ Fkl α
∂ Fim p −1 p α α −1 ∂ α α = F (t ) γ S0n j + Fim Fmn (t ) γ S0n j δn j − δ − ∂ Fkl mn ∂ Fkl n j α α
∂ γ α p p −1 α α −1 α δn j − = δik δml Fmn (t ) γ S0n j + Fim Fmn (t ) − S ∂ Fkl 0n j α α
p = δik Flnp (t )−1 δn j − γ α Sα0n j − Fim Fmn (t )−1 Rαkl Sα0n j , α
(B.5)
α
where
Rαkl =
∂ γ α . ∂ Fkl
• Calculation of Q:
Se (τ ) = Se trial (τ ) −
(B.6) α
γ α C α ( τ ) ,
Qi jkl =
∂Cα ∂ Siej ∂ Siejtrial α ∂ γ α = − Ci j − γ α i j = Di jkl − Ciαj Rαkl − γ α Jiαjkl , ∂ Fkl ∂ Fkl ∂ F ∂ F kl kl α α α α
Di jkl =
∂ Siejtrial , ∂ Fkl
(B.7)
(B.8)
where
Jiαjkl =
∂ Ciαj . ∂ Fkl
(B.9)
Let A = Ce trial (τ ) = Fe trial (τ )T Fe trial (τ ), thus,
Se trial = C : Ee trial = C :
Di jkl =
1 e trial 1 C − 1 = C : [A − 1], 2 2
1 ∂ Siejtrial ∂ 1 ∂ Amn 1 = Ci jmn (Amn − δmn ) = Ci jmn = Ci jmn Lmnkl , ∂ Fkl ∂ Fkl 2 2 ∂ Fkl 2
(B.10)
(B.11)
28
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
where
Lmnkl =
∂ Amn . ∂ Fkl
Recall
(B.12)
Cα = C : sym Ce trial Sα0 = C : sym[ASα0 ] =
1 T C : ASα0 + Sα0 A , 2
(B.13)
therefore,
∂ Ciαj 1 ∂ 1 = C A Sα + Sα0 pm A pn = Ci jmn Lmpkl Sα0 pn + Sα0 pm L pnkl . ∂ Fkl 2 ∂ Fkl i jmn mp 0 pn 2
Jiαjkl =
(B.14)
• Calculation of L: T
Fe trial (τ ) = F(τ )Fp (t )−1 ,
Ce trial (τ ) = Fe trial (τ ) Fe trial (τ ) = Fp (t )−T F(τ )T F(τ )Fp (t )−1 ,
(B.15)
∂ Ai j ∂ p −1 ∂ p p = Fmi (t ) Fnm Fnp Fppj (t )−1 = Fmi (t )−1 [F F ]F (t )−1 ∂ Fkl ∂ Fkl ∂ Fkl nm np p j p = Fmi (t )−1 δnk δml Fnp + Fnm δnk δ pl Fppj (t )−1
Li jkl =
p = Flip (t )−1 Fkm Fmp j (t )−1 + Fmi (t )−1 Fkm Flpj (t )−1
(B.16)
• Calculation of Rα :
γ α = γ˙ α (τ α , sαρ ) t = γ α (τ α , sαρ ), (B.17) where the slip rate γ˙ α is a function of the resolved stress τ α and the slip resistance sα . For the sake of simplicity, following ρ the assumption that the change in the slip resistance sα ρ during the time step is neglected (Balasubramanian, 1998; Dai, 1997), we have γ α = γ α (τ α ). For metallic materials, the elastic stretch is usually infinitesimal, therefore,
τ α = (Ce Se ) : Sα0 ≈ Se : Sα0 = τ α (Se ) ⇒ γ α = γ α (Se ),
(B.18)
and
Rαi j =
∂ γ α ∂ γ α ∂ Skle = = Bαkl Qkli j , ∂ Fi j ∂ Skle ∂ Fi j
Bαkl =
∂ γ α . ∂ Skle
(B.19)
Substituting the above equations into the calculation of Q, expressed in Eq. (B.8), we get:
Qi jkl = Di jkl − = Di jkl −
α
α
Ciαj Rαkl −
α
γ α Jiαjkl
Ciαj Bαmn Qmnkl −
or in tensorial form:
Q=D−
( Cα ⊗ Bα ) : Q −
α
α
α
γ α Jiαjkl ,
(B.20)
γ α J α .
(B.21)
Rearranging the terms in Eq. (B.21), we can solve for Q as follows:
Q=K
−1
:
D−
• Calculation of Bα :
Bα =
α
γ α J α
,
Ki jkl = δik δ jl +
α
Ciαj Bαkl .
(B.22)
∂ γ α ∂ γ α ∂τ α ∂ γ α ∂ e α ∂ γ α = = S : S = sym[Sα0 ]. ( ) 0 ∂ Se ∂τ α ∂ Se ∂τ α ∂ Se ∂τ α
(B.23)
Based on the derivations from Eq. (B.1) to Eq. (B.23), and recalling that and F = F(t + t ) are known, the computational procedure for calculating W (τ ), defined in Eq. (65), is summarized as follows: Fe (t),
p
p
p
p
1. Li jkl = Fli (t )−1 Fkm Fm j (t )−1 + Fmi (t )−1 Fkm Fl j (t )−1 , 2. Di jkl = 12 Ci jmn Lmnkl , 3. Jiαjkl = 12 Ci jmn [Lmpkl Sα + Sα L ], 0 pn 0 pm pnkl ∂ γ α
4. Bα = ∂τ α sym[Sα ], ij 0i j 5. Ki jkl = δik δ jl + α Ciαj Bα , kl α α 6. Qi jkl = Ki−1 ( D − mnkl α γ Jmnkl ), jmn αQ , 7. Rα = B ij kl kli j p p 8. Ti jkl = δik Fln (t )−1 (δn j − α γ α Sα ) − Fim Fmn (t )−1 α Rαkl Sα0n j , 0n j 9. Wi jkl =
1 [T Se F e detFe imkl mn jn
e Q e e e e e e e −1 + Fim mnkl Fjn + Fim Smn T jnkl − Fim Smn Fjn Fqp T pqkl ].
Fp (t),
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
29
References Acharya, A., Beaudoin, A., 20 0 0. Grain-size effect in viscoplastic polycrystals at moderate strains. J. Mech. Phys. Solids 48 (10), 2213–2230. Addessio, F.L., Luscher, D.J., Cawkwell, M.J., Ramos, K.J., 2017. A single-crystal model for the high-strain rate deformation of cyclotrimethylene trinitramine including phase transformations and plastic slip. J. Appl. Phys. 121 (18), 185902. Agrawal, V., Dayal, K., 2015. A dynamic phase-field model for structural transformations and twinning: regularized interfaces with transparent prescription of complex kinetics and nucleation. part i: formulation and one-dimensional characterization. J. Mech. Phys. Solids 85, 270–290. Agrawal, V., Dayal, K., 2015. A dynamic phase-field model for structural transformations and twinning: regularized interfaces with transparent prescription of complex kinetics and nucleation. part II: two-dimensional characterization and boundary kinetics. J. Mech. Phys. Solids 85, 291–307. Anand, L., 2004. Single-crystal elasto-viscoplasticity: application to texture evolution in polycrystalline metals at large strains. Comput. Methods Appl. Mech. Eng. 193 (48), 5359–5383. Ardeljan, M., Beyerlein, I., McWilliams, B., Knezevic, M., 2016. Strain rate and temperature sensitive multi-level crystal plasticity model for large plastic deformation behavior: application to AZ31 magnesium alloy. Int. J. Plast. 83, 90–109. Ardeljan, M., McCabe, R., Beyerlein, I., Knezevic, M., 2015. Explicit incorporation of deformation twins into crystal plasticity finite element models. Comput. Methods Appl. Mech. Eng. 295, 396–413. Asaro, R., Needleman, A., 1985. Texture development and strain hardening in rate dependent polycrystals. Acta Metall. 33, 923–953. Asaro, R., Rice, J., 1977. Strain localization in ductile single crystals. J. Mech. Phys. Solids 25 (5), 309–338. Asaro, R.J., 1983. Micromechanics of crystals and polycrystals. In: Advances in Applied Mechanics, 23. Elsevier, pp. 1–115. Balasubramanian, S., 1998. Polycrystalline Plasticity: Application to Deformation Processing of Lightweight Metals. Massachusetts Institute of Technology. Belytschko, T., Fish, J., Engelmann, B.E., 1988. A finite element with embedded localization zones. Comput. Methods Appl. Mech. Eng. 70, 59–89. Belytschko, T., Liu, W.K., Moran, B., Elkhodary, K., 2013. Nonlinear finite elements for continua and structures.. John wiley & sons. Beyerlein, I., Capolungo, L., Marshall, P., McCabe, R., Tomé, C., 2010. Statistical analysis of deformation twinning in magnesium. Phil. Mag. 90, 2161–2190. Beyerlein, I., Tomé, C., 2010. A probabilistic twin nucleation model for HCP polycrystalline metals. Proc. Royal Soc. A 466, 2517–2544. Beyerlein, I.J., Tomé, C.N., 2008. A dislocation-based constitutive law for pure Zr including temperature effects. Int. J. Plast. 24, 867–895. Bronkhorst, C., Hansen, B., Cerreta, E., Bingert, J., 2007. Modeling the microstructural evolution of metallic polycrystalline materials under localization conditions. J. Mech. Phys. Solids 55 (11), 2351–2383. Busso, E.P., 1990. Cyclic deformation of monocrystalline nickel aluminide and high temperature coatings. Massachusetts Institute of Technology. Capolungo, L., Marshall, P., McCabe, R., Beyerlein, I., Tomé, C., 2009. Nucleation and growth of twins in Zr: a statistical study. Acta Mat. 57, 6047–6056. Cheng, J., Ghosh, S., 2017. Crystal plasticity finite element modeling of discrete twin evolution in polycrystalline magnesium. J. Mech. Phys. Solids 99, 512–538. Cheng, J., Shen, J., Mishra, R.K., 2018. Discrete twin evolution in Mg alloys using a novel crystal plasticity finite element model. Acta Mat. 149, 142–153. Chester, S.A., Bernier, J.V., Barton, N.R., Balogh, L., Clausen, B., Edmiston, J.K., 2016. Direct numerical simulation of deformation twinning in polycrystals. Acta Mat. 120, 348–363. Clayton, J.D., 2011. Nonlinear Mechanics of Crystals. Springer. Clayton, J.D., Knap, J., 2011. A phase field model of deformation twinning: nonlinear theory and numerical simulations. Physica D 240, 841–858. Clayton, J.D., Knap, J., 2016. Phase field modeling and simulation of coupled fracture and twinning in single crystals and polycrystals. Comput. Methods Appl. Mech. Eng. 312, 447–467. Dai, H., 1997. Geometrically-necessary dislocation density in continuum plasticity theory, FEM implementation and applications. Massachusetts Institute of Technology. Dequiedt, J.L., Denoual, C., Madec, R., 2015. Heterogeneous deformation in ductile FCC single crystals in biaxial stretching: the influence of slip system interactions. J. Mech. Phys. Solids 83, 301–318. Feng, B., Bronkhorst, C., Addessio, F., Morrow, B., Cerreta, E., Lookman, T., Lebensohn, R., Low, T., 2018. Coupled elasticity, plastic slip, and twinning in single crystal titanium loaded by split-hopkinson pressure bar. J. Mech. Phys. Solids 119, 274–297. Fish, J., Belytschko, T., 1988. Elements with embedded localization zones for large deformation problems. Comp. Struct. 30, 247–256. Ghosh, S., Cheng, J., 2018. Adaptive multi-time-domain subcycling for crystal plasticity FE modeling of discrete twin evolution. Comput. Mech. 61, 33–54. Gray III, G.T., 1997. Influence of strain rate and temperature on the structure-property behavior of high-purity titanium. J. Phys. IV France 7. C3–423–C3–428 Gurao, N.P., Kapoor, R., Suwas, S., 2011. Deformation behavior of commercially pure titanium at extreme strain rates. Acta Mat. 59 (9), 3431–3446. Gurtin, M.E., Fried, E., Anand, L., 2010. The Mechanics and Thermodynamics of Continua. Cambridge University Press. Hansen, B.L., Bronkhorst, C.A., Ortiz, M., 2010. Dislocation subgrain structures and modeling the plastic hardening of metallic single crystals. Modelling Simul. Mater. Sci. Eng. 18, 055001. Hill, R., Rice, J., 1972. Constitutive analysis of elastic-plastic crystals at arbitrary strain. J. Mech. Phys. Solids 20 (6), 401–413. Hughes, T.J.R., 20 0 0. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover Publications. Jin, T., Mourad, H.M., Bronkhorst, C.A., 2019. A comparative study of shear band tracking strategies in three-dimensional finite elements with embedded weak discontinuities. Finite Elem. Anal. Des. 155, 11–31. Jin, T., Mourad, H.M., Bronkhorst, C.A., Livescu, V., 2018. Finite element formulation with embedded weak discontinuities for strain localization under dynamic conditions. Comput. Mech. 61, 3–18. Jin, T., Mourad, H.M., Bronkhorst, C.A., Livescu, V., Zhang, X., Linder, C., Regueiro, R.A., 2019. Three-dimensional explicit finite element formulation for shear localization with global tracking of embedded weak discontinuities. Comput. Methods Appl. Mech. Eng. 353, 416–447. Kalidindi, S., Bronkhorst, C., Anand, L., 1992. Crystallographic texture evolution in bulk deformation processing of FCC metals. J. Mech. Phys. Solids 40 (3), 537–569. Kalidindi, S.R., 1998. Incorporation of deformation twinning in crystal plasticity models. J. Mech. Phys. Solids 46 (2), 267–290. Knezevic, M., Lebensohn, R., Cazacu, O., Revil-Baudard, B., Proust, G., Vogel, S.C., Nixon, M., 2013. Modeling bending of α -titanium with embedded polycrystal plasticity in implicit finite elements. Mater. Sci. Eng. A 564, 116–126. Kocks, U., 1976. Laws for work-hardening and low-temperature creep. J. Eng. Mater. Technol. 98 (1), 76–85. Kothari, M., Anand, L., 1998. Elasto-viscoplastic constitutive equations for polycrystalline metals: application to tantalum. J. Mech. Phys. Solids 46, 51–83. Kröner, E., Teodosiu, C., 1974. Lattice defect approach to plasticity and viscoplasticity. In: Sawczuk, A. (Ed.), Problems of plasticity. Noordhoff International Publishing, Leiden, pp. 45–82. Kubin, L.P., 2013. Dislocations, mesoscale simulations and plastic flow. Oxford University Press. Kumar, M.A., Beyerlein, I.J., McCabe, R.J., Tomé, C.N., 2016. Grain neighbor effects on twin transmission in hexagonal close-packed materials. Nature Commun. 7, 13826. Kumar, M.A., Beyerlein, I.J., Tomé, C.N., 2016. Effect of local stress fields on twin characteristics in HCP metals. Acta Mat. 116, 143–154. Kumar, M.A., Capolungo, L., McCabe, R.J., Tomé, C.N., 2019. Characterizing the role of adjoining twins at grain boundaries in hexagonal close packed materials. Sci. Rep. 9, 3846. Kumar, M.A., Kanjarla, A.K., Niezgoda, S.R., Lebensohn, R.A., Tomé, C.N., 2015. Numerical study of the stress state of a deformation twin in magnesium. Acta Mat. 84, 349–358. Leclercq, S., Nguy, C., Bensussan, P., 1989. Microstructure transformations in pure polycrystalline α -titanium under static and dynamic loading. Mechanical Properties of Materials at High Rates of Strain. Institute of Physics, Bristol. Lee, E., 1969. Elastic-plastic deformation at finite strains. J. Appl. Mech. 36 (1), 1–6. Lee, E.H., Liu, D.T., 1967. Finite-strain elastic–plastic theory with application to plane-wave analysis. J. Appl. Phys. 38 (1), 19–27.
30
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Levitas, V.I., Preston, D.L., 2005. Thermomechanical lattice instability and phase field theory of martensitic phase transformations, twinning and dislocations at large strains. Phys. Lett. A 343, 32–39. Liu, C., Shanthraj, P., Diehl, M., Roters, F., Dong, S., Dong, J., Ding, W., Raabe, D., 2018. An integrated crystal plasticity-phase field model for spatially resolved twin nucleation, propagation, and growth in hexagonal materials. Int. J. Plast. 106, 203–227. Livescu, V., Beyerlein, I.J., Bronkhorst, C.A., Dippo, O.F., Ndefru, B.G., Capolungo, L., Mourad, H.M., 2019. Microstructure insensitive twinning: a statistical analysis of incipient twins in high-purity titanium. Materialia 6, 100303. Luan, Q.M., Britton, T.B., Jun, T.S., 2018. Strain rate sensitivity in commercial pure titanium: the competition between slip and deformation twinning. Mater. Sci. Eng. A 734, 385–397. Mourad, H.M., Bronkhorst, C.A., Livescu, V., Plohr, J.N., Cerreta, E.K., 2017. Modeling and simulation framework for dynamic strain localization in elasto-viscoplastic metallic materials subject to large deformations. Int. J. Plast. 88, 1–26. Niezgoda, S.R., Kanjarla, A.K., Beyerlein, I.J., Tomé, C.N., 2014. Stochastic modeling of twin nucleation in polycrystals: an application in hexagonal close– packed metals. Int. J. Plast. 56, 119–138. Peirce, D., Asaro, R., Needleman, A., 1982. An analysis of nonuniform and localized deformation in ductile single crystals. Acta Metall. 30 (6), 1087–1119. Rapperport, E., Hartley, C., 1960. Deformation modes of zirconium at 77 K, 300 K, 575 K, and 1075 K. Trans. Am. Inst. Min. Metall. Eng. 218, 869. Rice, J., 1971. Inelastic constitutive relations for solids: an internal-variable theory and its application to metal plasticity. J. Mech. Phys. Solids 19 (6), 433–455. Salem, A.A., Kalidindi, S.R., Doherty, R.D., 2003. Strain hardening of titanium: role of deformation twinning. Acta Mat. 51 (14), 4225–4237. Salem, A.A., Kalidindi, S.R., Doherty, R.D., Semiatin, S.L., 2006. Strain hardening due to deformation twinning in α -titanium: mechanisms. Metall. Mater. Trans. A 37 (1), 259–268. Simo, J.C., Hughes, T.J.R., 1986. On the variational foundations of assumed strain methods. J. Appl. Mech. 53, 51–54. Song, S., III, G.G., 1995. Influence of temperature and strain rate on slip and twinning behavior of Zr. Metall. Mater. Trans. A 26 (10), 2665–2675. Tjahjanto, D.D., Turteltaub, S., Suiker, A.S.J., 2008. Crystallographically based model for transformation-induced plasticity in multiphase carbon steels. Continuum Mech. Thermodyn. 19 (7), 399–422. Van Houtte, P., 1978. Simulation of the rolling and shear texture of brass by the Taylor theory adapted for mechanical twinning. Acta Metall. 26 (4), 591–604. Wang, J., Hirth, J., Tomé, C., 2009. (1¯ 012 ) twinning nucleation mechanisms in hexagonal-close-packed crystals. Acta Mat. 57, 5521–5530. Wang, J., Hoagland, R., Hirth, J., Capolungo, L., Beyerlein, I., Tomé, C., 2009. Nucleation of a (1¯ 012 ) twin in hexagonal close-packed crystals. Scripta Mater. 61, 903–906. Wang, J., Beyerlein, I.J, Tomé, C., 2010. An atomic and probabilistic perspective on twin nucleation in Mg. Scripta Mater. 63, 741–746. Wang, J., Yadav, S., Hirth, J., Tomé, C., 2013. Pure-shuffle nucleation of deformation twins in hexagonal-close-packed metals. Mater. Res. Lett. 1, 126–132.
Contents lists available at ScienceDirect
Journal of the Mechanics and Physics of Solids journal homepage: www.elsevier.com/locate/jmps
A single crystal plasticity finite element formulation with embedded deformation twins Tao Jin a, Hashem M. Mourad a,∗, Curt A. Bronkhorst b, Irene J. Beyerlein c a
Theoretical Division, T-3, Los Alamos National Laboratory, Los Alamos, NM 87545, USA Department of Engineering Physics, University of Wisconsin, Madison, WI, 53706, USA c Department of Mechanical Engineering, Materials Department, University of California, Santa Barbara, CA, 93106, USA b
a r t i c l e
i n f o
Article history: Received 31 May 2019 Revised 23 August 2019 Accepted 10 September 2019 Available online 14 September 2019 Keywords: Deformation twinning Embedded weak discontinuity Single-crystal elasto-viscoplasticity
a b s t r a c t Deformation twinning is an important plastic deformation mechanism in some polycrystalline metals such as titanium and magnesium. In this paper, we present a novel crystal plasticity finite element framework that accounts for deformation twinning explicitly, in addition to crystallographic slip. Within this computational framework, deformation twins are treated as weak discontinuities embedded within individual finite elements, such that a jump in the velocity gradient field is introduced (via the discretized gradient operator) between the twinned and untwinned crystalline regions, taking into account compatibility and traction continuity conditions at the interface between these two regions. The deformation gradient is multiplicatively split into elastic and plastic parts in the untwinned region, as is customary in finite-deformation crystal plasticity formulations. A different multiplicative decomposition of the deformation gradient into elastic, plastic (slip), and twinning parts is adopted in the twinned region, allowing deformation twinning to be accounted for as an additional mode of plastic deformation. A stochastic model is used to predict twin nucleation at grain boundaries, and the evolution of the length and thickness of the twinned region under subsequent deformation is taken into account. The linearization of the single-crystal plasticity model, required in order to enforce traction continuity at the interface between the twinned and untwinned regions using a Newton iterative scheme, is described in detail. Simulations of the initiation and propagation of {101¯ 2} tensile twins in hexagonal close packed (HCP) titanium are presented to demonstrate the capabilities of the proposed computational framework. © 2019 Elsevier Ltd. All rights reserved.
1. Introduction The plastic deformation of metallic materials is in general mediated by the motion of dislocations, and the creation and propagation of deformation twins under mechanical loading. In order to accommodate imposed deformation, a metallic material will naturally select combinations of dislocation motion and deformation twinning. The propensity for the natural combination of the two plastic deformation mechanisms depends on the material’s crystallographic structure as well as the boundary conditions imposed upon the material to produce a given deformation field. High-symmetry materials such as copper (FCC) and tantalum (BCC) have several crystallographically close-packed atomic planes whose orientations are ∗
Corresponding author. E-mail addresses: [email protected] (H.M. Mourad), [email protected] (C.A. Bronkhorst).
https://doi.org/10.1016/j.jmps.2019.103723 0022-5096/© 2019 Elsevier Ltd. All rights reserved.
2
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
evenly spaced in three-dimensional (3D) space, allowing them to accommodate complex deformation fields via plastic slip (Clayton, 2011). Generally, the initial mobile dislocation content in such high-symmetry materials is sufficiently high, so that the resolved shear stress necessary to initiate dislocation motion is small relative to that required to nucleate deformation twins (Kubin, 2013). In addition, the Burgers vector magnitude is uniform for all systems in these crystal structures, and therefore, there is not a significant bias in resistance to dislocation motion (Bronkhorst et al., 2007; Kothari and Anand, 1998). Deformation twinning is still observed in these materials under certain loading conditions, however, generally, it is not a prominent feature. Materials such as titanium (HCP) are of low-symmetry crystallographic structure. Members of this class of materials possess a number of close-packed atomic planes, but their atomic unit cell is hexagonally cylindrical, and therefore the available slip planes/directions are not uniformly spaced in 3D. When loaded in certain directions, the unit cell does not have an adequate number of slip systems available to accommodate plastic deformation (Ardeljan et al., 2016; Beyerlein and Tomé, 2008; Gray III, 1997; Knezevic et al., 2013; Rapperport and Hartley, 1960; Salem et al., 20 03; 20 06; Song and Gray, 1995). In addition, the Burgers vector magnitude is significantly different for each of the three categories of slip systems present in HCP materials, and therefore, the resolved shear stress needed to give rise to dislocation motion also varies correspondingly (a slip system with a larger Burgers vector magnitude generally presents a greater resistance to dislocation glide; e.g. see Liu et al., 2018). For HCP materials, some estimates put these differences at several times the minimum value. The probability that deformation twinning will be required to accommodate an imposed deformation is then greater for HCP materials, and this is frequently observed experimentally (Beyerlein et al., 2010; Capolungo et al., 2009; Luan et al., 2018). Within the field of computational plasticity, deformation twinning offers unique challenges in comparison with dislocation slip. Although the physics of both mechanisms are discrete by nature, the characteristic size of deformation twins is much greater than that of dislocations (Gray, 1997; Gurao et al., 2011; Leclercq et al., 1989). Given that computational tools such as the finite element method rely upon spatial discretization of a body into individual computational elements, this disparity in characteristic size presents a challenge. Accounting for both mechanisms in full-field computational approaches such as the crystal plasticity finite element (CPFE) method therefore requires a novel approach. Although still a vigorous research area, it has become standard to represent the mechanics of dislocation motion using continuum single crystal plasticity theories (Clayton, 2011; Gurtin et al., 2010; Kubin, 2013). The characteristic length scale of dislocations still allows for an approximate homogenized representation within CPFE approaches. It still remains a challenge to adequately represent the complex interactions between dislocations and the subcell structures developed with deformation (e.g. Dequiedt et al., 2015; Hansen et al., 2010). With the characteristic length scale of deformation twins being comparable to the grain size, it is difficult to be content with a homogenized representation at the same level as for dislocations. Within the CPFE context, there is a distinct need to provide a suitable computational framework for deformation twinning, which also accounts for the evolving interfaces created. The goal of the present work is to provide this kind of differentiated representation of dislocations and deformation twins, within the context of CPFE. The process of deformation twinning is believed to involve three steps: nucleation, propagation, and growth (Beyerlein and Tomé, 2010; Cheng and Ghosh, 2017; Livescu et al., 2019; Niezgoda et al., 2014). For relatively high-purity materials, grain boundary (GB)–dislocation interaction is hypothesized to be one possible mechanism of twin nucleation (Wang et al., 20 09a; 20 09b; 2010; 2013). In this mechanism, dislocations impinging upon a GB1 interact with the latter to produce a number of partial dislocations at the impingement site (Niezgoda et al., 2014). This localized partial dislocation nucleus must be of a certain size or energy level in order for a twin to be formed within the target single crystal (Liu et al., 2018; Niezgoda et al., 2014; Wang et al., 20 09a; 20 09b; 2010; 2013). This initiation process relies upon the storage of sufficient energy within the nucleus at the GB, to allow it to propagate into the neighboring grain. In simplest form, initiation leads to the growth stage as the boundaries of the twinned region migrate in the width direction. The corresponding transformation has both shear deformation and crystallographic rotation kinematic components across the moving interface surface (Ardeljan et al., 2015; Cheng and Ghosh, 2017; Cheng et al., 2018). There have been several recent efforts to extend computational treatments of deformation twinning beyond homogenized representations. This is necessitated not only by the need to improve modeling of plasticity in general, but also by the need to integrate computational and experimental work in order to explore the science of deformation twinning more effectively. Within the domain of computational meso-mechanics for explicit representation of twinning, there has been important recent development work. Phase field techniques have recently begun to be used to examine the local evolution in structure with the inclusion of deformation twinning. The coupled processes of martensitic phase transformation together with both dislocation motion and twinning were presented by Levitas and Preston (2005) within a formal thermodynamic framework targeting NiAl, BN, C materials. The theory was formulated for large deformation; however, limited computational results were presented. A different model that represents phase transformation and deformation twinning was presented in a series of two papers (Agrawal and Dayal, 2015a; 2015b), representing a progression from one to two dimensions. This work focused on the advanced treatment of evolving interfaces associated with both transformation and twinning, represented implicitly by a continuous order parameter field. Two other publications (Clayton and Knap, 2011; 2016) firstly presented an advanced formulation for anisotropic single crystal materials along with anisotropic interface representation for twin boundaries. The
1 It is important to keep in mind that each GB is atomistically non-uniform and that dislocation impingement upon a GB surface will be discrete and non-uniform (Wang et al., 2009a, 2009b, 2010, 2013), although there is still much to be learned about the specifics.
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
3
numerical implementation was described and two-dimensional proof-of-principle calculations were also presented. In the following paper, the authors progressed to 3D polycrystal systems, while including the physics of brittle damage. Coupled twinning and damage calculations were presented to demonstrate the importance and potential of advanced meso-scale computational capability. HCP materials were the target of a recent work (Liu et al., 2018) where high-fidelity phase field representation, coupled with classical models for dislocation slip, were applied to polycrystal systems of magnesium. That work paid particular attention to nucleation at GBs, which was examined in the context of the nucleation, propagation and growth progression of interest here. HCP materials were also the focus of a series of papers (Kumar et al., 2016a; 2016b; 2019; 2015) where the spectral FFT solver was used to examine the behavior of tri-crystal configurations. Traditional finite element tools were employed by three groups to provide an explicit spatial representation of twins in 3D (Ardeljan et al., 2015; Cheng and Ghosh, 2017; Cheng et al., 2018; Chester et al., 2016; Ghosh and Cheng, 2018). This achieved significant advancements in relation to the traditional CPFE method and articulated the resolution challenges for this evolving length scale problem. The computational framework proposed here is an extension of one that has been developed previously for the representation of shear banding (Jin et al., 2019a; 2018; 2019b; Mourad et al., 2017). Recognizing that the topology of shear bands and deformation twins is similar, this computational framework is extended here to the kinematics and physics of large deformation within metallic single crystals. The current work differs from the aforementioned ones in that we propose an alternative weak-form embedding technique to explicitly represent individual deformation twins and their progression. The remainder of this paper is organized as follows. In Section 2.1, we describe the kinematics associated with the twinned and untwinned regions in the crystal and the multiplicative decomposition of the deformation gradient in each region. In Section 2.2, we discuss in detail the embedded weak discontinuity approach, which is adopted to treat deformation twinning as a heterogeneous deformation mode. In Section 2.3, we present the single-crystal hyperelastic–viscoplasticity model, and we introduce the probability-based twin nucleation model and the evolution equations of the twin growth in Section 2.4. In Section 3, we present several numerical examples to simulate the initiation and propagation of the {101¯ 2} tensile twin system in hexagonal close packed (HCP) titanium. Finally, a summary of the present paper and possible future work directions are discussed in Section 4. 2. Computational framework 2.1. Kinematics We consider a crystalline body with the reference configuration 0 undergoing the deformation ϕ over a time interval of interest T . At time t ∈ T , every particle X ∈ 0 is mapped to a corresponding point x in the deformed configuration , such that
x = ϕ ( X , t ).
(1)
The deformation gradient tensor is defined as follows:
F :=
∂ϕ = FiJ ei ⊗ E J , ∂X
(2)
where {EI } and {ei } denote the Cartesian basis vectors in the reference and current configurations, respectively. Using superposed dots to designate material time derivatives, the velocity gradient is expressed as
L :=
∂ v(x, t ) ˙ −1 = FF , ∂x
(3)
where v := u˙ is the velocity vector, defined in terms of the displacement u := x − X . Assuming that the motion of dislocations takes place on N prescribed slip systems in the crystalline lattice, the α -th α slip system (α = 1, 2, . . . , N ) is described by a pair of orthogonal unit directional vectors {t α 0 , s0 } in the reference configuα represents the slip direction. The evolution of the plastic ration, where t α represents the normal of the slip plane, and s 0 0 deformation gradient Fp is given by
F˙ p = Lp Fp ,
Lp =
α
γ˙ α Sα0 =
α
γ˙ α sα0 ⊗ t α0 ,
(4)
α where γ˙ α is the plastic slip rate on the α -th slip system. The slip plane normal t α 0 , the slip direction s0 , and the Schmid α ⊗ t α , are pushed forward according to the following relationships, tensor Sα := s 0 0 0
t α = Fe−T t α0 ,
sα = Fe sα0 ,
Sα = Fe Sα0 Fe
−1
= sα ⊗ t α ,
(5)
where tα , sα and Sα are the corresponding quantities in the deformed configuration, and Fe := FFp −1 is the elastic deformation gradient obtained via multiplicative decomposition of F, which leads to the following additive decomposition of the velocity gradient:
L = Le +
α
γ˙ α Sα , Le = F˙ e Fe −1 .
(6)
4
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Fig. 1. Multiplicative decomposition of the deformation gradient FU for the untwinned region into the elastic part FU,e and the plastic part FU,p .
Fig. 2. Multiplicative decomposition of the deformation gradient tensor FZ for the twinned region into the elastic part FZ,e , the plastic slip part FZ,p , and the twinning part FZ,tw .
After deformation twinning is initiated, plastic deformation is accommodated by both plastic slip and twinning. In this paper, we treat deformation twinning as a heterogeneous deformation mode, and represent it using the embedded weakdiscontinuity approach detailed in Section 2.2. We use (•)U and (•)Z to represent quantities associated with the untwinned and twinned crystal regions, respectively. Thus, in the untwinned region where plastic deformation is accommodated by dislocation glide on the prescribed slip systems, the deformation gradient FU is decomposed, as illustrated in Fig. 1, in accordance with the classical multiplicative decomposition (Kröner and Teodosiu, 1974; Lee, 1969; Lee and Liu, 1967):
FU = FU,e FU,p ,
(7)
where FU,p characterizes the local plastic deformation due to dislocation motion and FU,e represents the elastic stretching and rotation of the crystal lattice. The plastic part of the deformation gradient, FU,p , evolves according to Eq. (4); i.e.
F˙ U,p = LU,p FU,p ,
LU,p =
α
γ˙ α Sα0 =
α
γ˙ α sα0 ⊗ t α0 .
(8)
In the twinned region, deformation twinning is included as an additional mode of plastic deformation, as illustrated in Fig. 2. As a result, the multiplicative decomposition of the deformation gradient takes the following form (Addessio et al., 2017; Feng et al., 2018; Tjahjanto et al., 2008):
FZ = FZ,e FZ,p FZ,tw .
(9) FZ,tw ,
The twinning part of the deformation gradient, which maps the reference configuration onto the twinning configuration (the first relaxed intermediate configuration shown in Fig. 2), is given by η
FZ,tw = 1 + γtw mη ⊗ nη ,
(10) η
where 1 is the second-order identity tensor, and γtw , mη and nη are, respectively, the magnitude of the twinning shear
strain, the shearing direction, and the unit normal to the habit plane, characteristic of the (η-th) twin system under consideration. Furthermore, in Eq. (9), the elastic stretching and rotation of the crystalline lattice inside the twin is represented by FZ,e . Plastic deformation due to dislocation motion inside the twinned region is characterized by FZ,p , which evolves according to
F˙ Z,p = LZ,p FZ,p ,
LZ,p =
α
α
γ˙ α Sˆ 0 =
α
α
γ˙ α sˆ α0 ⊗ tˆ 0 .
(11) α
α
Here, LZ,p represents the plastic velocity gradient in the twinned region. The slip plane normal tˆ 0 , the slip direction sˆ 0 and α α ˆ 0 := sˆ α ˆ the Schmid tensor S 0 ⊗ t 0 , corresponding to the α -th slip system in the twinned region of the crystal, are related to
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
5
their counterparts in the untwinned region via α
α
tˆ 0 = R t α0 ,
sˆ 0 = R sα0 ,
α
ˆ 0 = R Sα0 RT , S
(12)
where R is an orthogonal matrix representing the rotation between the lattice orientations in the untwinned and twinned regions, given by (Kalidindi, 1998; Van Houtte, 1978):
R = 2 n 0 ⊗ n 0 − 1,
(13)
in terms of the unit vector n0 , normal to the twinning plane in the reference configuration. It is important to note that the tensorial quantities used to describe the deformation in the untwinned and twinned regions are not completely independent as they need to satisfy compatibility conditions at the interface formed between these two regions. In addition, the stresses in the crystal must satisfy traction continuity (equilibrium) conditions at the parent– twin interface. In Section 2.2, we discuss these constraints in more detail and describe the methods used to enforce them within the present computational framework, where deformation twins are represented as embedded weak discontinuities. 2.2. Embedded weak discontinuity finite element formulation Representation of deformation twinning in a computational setting presents unique challenges compared to dislocation slip, mainly owing to the large disparity in length scales associated with these two plastic deformation mechanisms. While the length scale associated with dislocations allows for an approximate homogenized representation in CPFE simulations, it is difficult to be content with a volume fraction-based, pseudo-slip representation of twinning, in view of the fact that its characteristic length scale is comparable to the grain size. One limitation of such homogenized treatments of twinning, for instance, is related to their inability to resolve the local stress state at the twin–parent grain interface, driving the growth of an existing twin. The need for a computational approach which can address such limitations has led to the development of CPFE approaches where twins are resolved explicitly using conventional finite element methods coupled with single crystal plasticity models (Ardeljan et al., 2015). Such an approach, however, requires remeshing of every grain containing a twin, immediately after the twin is formed and at every subsequent time step while the twin is growing, which can render CPFE simulations prohibitively expensive, especially since a highly refined mesh is often needed to discretize twinned regions due to their typical high aspect ratio. To address the aforementioned computational challenges, we adopt an embedded weak-discontinuity finite element formulation, based on a generalization of the Hu–Washizu variational principle to dynamic problems, where the velocity, velocity gradient, and stress are treated as three independent fields. This allows twinning to be treated as a heterogeneous deformation mode, represented within an individual finite element as a discontinuity in the velocity gradient field (i.e. a weak discontinuity). In the present section, we show how the assumed velocity gradient, which appears in the Hu–Washizu variational statement of the problem, is designed to admit such weak discontinuities, such that twins can be embedded within the finite element grid (mesh) instead of being resolved explicitly. This leads to a sub-grid finite element formulation, capable of representing features of the solution (in this case, weak discontinuities representing deformation twinning) that are smaller than the size of the finite element grid, which obviates the need for highly refined meshes and expensive remeshing operations. Given that twinning tends to be more pronounced in dynamic problems involving high-strain rate loading conditions, and since the twin propagation velocity is comparable to the shear wave speed in the material, requiring extremely small time steps in order to resolve the relevant time scales, the formulation includes inertial effects and adopts an explicit time stepping scheme. 2.2.1. Variational form Introducing the assumed velocity gradient L and the assumed stress S, and treating these as independent field variables along with the velocity v, the three-field variational form is expressed as follows,
δ =
δv · (ρ v˙ − ∇ · S ) d +
+
δ S : ( L − L ) d +
h
δ L : ( σ − S ) d
δv · (S · n¯ − h¯ ) d = 0,
(14)
revealing the corresponding Euler–Lagrange equations,
ρ v˙ − ∇ · S = 0 S = σ ( L, q ) L=L S · n¯ = h¯
(balance of linear momentum), (constitutive relation), (compatibility condition), (Neumann boundary condition).
(15)
Here, h¯ is the applied traction acting upon the subset h of the domain’s boundary, which has an outward unit normal vector denoted by n¯ . Recall that L = ∇ v is the velocity gradient, as defined in Eq. (3) and note that σ = σ (L, q ) is the stress obtained through the constitutive relations, where q denotes internal state variables. It is also noted that Eq. (14) is the same updated Lagrangian weak form used for the treatment of shear banding in previous works (Jin et al., 2019a; 2018; 2019b; Mourad et al., 2017). While it is possible to develop an alternative total Lagrangian weak form, this additional development
6
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
effort is not undertaken herein since the two Lagrangian formulations are essentially identical in terms of accuracy and computational performance (see Belytschko et al., 2013, pp. 147–148). 2.2.2. Interpolation scheme The deformed configuration is partitioned into ne disjoint element domains, e , with boundaries e , and the following interpolation scheme for the fields v, L, L, and S, is adopted:
vi = NiA d˙ A ,
(16)
Li j = B¯ i jA d˙ A ,
(17)
˜ i jA d˙ A , Li j = B
(18)
Si j = Si jA sA ,
(19)
where d˙ A denotes a velocity degree of freedom,2 NiA is the corresponding nodal shape function, sA is a stress degree of ¯ operator in freedom, Si jA is a stress interpolation function, and summation is implied over repeated indices. Using the B the discretization of L is equivalent to selective reduced integration of volumetric and shear strain terms (see Mourad et al., 2017), and thus alleviates volumetric locking arising in the presence of well-developed (isochoric) plastic flow, and mitigates parasitic shear effects manifested under bending-dominated conditions. It is noted that the assumed velocity gradient L is ˜ , described in detail below. discretized using a different operator, B 2.2.3. Finite element equations Substituting Eqs. (16)–(19) into (14) and writing every integral over the global domain (or the subset h of its boundary) as a sum of integrals over each element e (or subset he of the element’s boundary), we obtain
δ d˙ A
MAB d¨ B + fAint − fAext + δ sA HAB d˙ B = 0,
(20)
for all arbitrary variations δ d˙ A and δ sA , i.e., MAB d¨ B + fAint = fAext ,
(21)
HAB d˙ B = 0,
(22)
where MAB =
ne
e
e=1
HAB =
ne
fAint =
e
ne e=1
fAext =
e=1
e
ne
he
e=1
ρ NiA NiB d,
(23)
¯ i jB − B ˜ i jB d, Si jA B
(24)
˜ i jA σi j d + HBA sB , B
(25)
NiA h¯ i d .
(26)
As shown by Simo and Hughes (1986), ensuring that Eq. (22) holds for any admissible velocity field requires (L − L ) to be orthogonal to the admissible stress field S, i.e. H = 0. This orthogonality condition is expressed as follows:
e
¯ −B ˜ d = 0, B
(27)
assuming a piecewise constant stress field. As a result, Eq. (25) is simplified, and the internal force is given by fAint =
ne e=1
e
˜ i jA σi j d. B
(28)
It should be noted here that σ ij is the stress obtained through the constitutive equations as a function of the assumed velocity gradient L. 2 Here we adopt the notation of Belytschko et al. (1988); Fish and Belytschko (1988). This notation reflects our computer implementation where degrees of freedom are stored in one-dimensional arrays. For example, individual components of nodal velocity vectors are stored in the one-dimensional array d˙ , and the index of this array is denoted by a capital letter, A or B, and encapsulates both the node number and the component number (i.e. direction in space) corresponding to each nodal degree of freedom.
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
7
Fig. 3. In the present approach, a deformation twin is modeled as a weak discontinuity embedded within individual finite elements; i.e. a strain jump is introduced across the interface between the embedded twin RZ and the surrounding parent grain (untwinned region) RU .
2.2.4. Assumed velocity gradient field In the embedded weak discontinuity approach, the velocity (and displacement) fields remain continuous, while discontinuities are introduced into the velocity gradient (strain rate) field to describe strain jumps across the interface between the embedded twin and the surrounding material (untwinned region), as shown schematically in Fig. 3. The present embedded discontinuity approach is classified as a sub-grid method as it allows a twinned region of width b ≤ h to be embedded within a finite element with characteristic size h. The jump in the velocity gradient field assumes the following mathematical form (Fish and Belytschko, 1988):
[[L]] = g m ⊗ n = g Ti j ei ⊗ e j ,
(29)
where g is the magnitude of the jump, n is the unit vector normal to the twin plane, and m is a unit vector such that T = m ⊗ n represents the discontinuous deformation mode. Noting that L: T represents the projection of the velocity gradient onto the discontinuous mode, the assumed velocity gradient is written in component form as follows:
Li j = Li j + gTi j =
δik δ jl + α (x )Ti j Tkl Lkl ,
(30)
αZ
where δ ij is the Kronecker delta, and α (x) is a scalar parameter such that α (x ) = in the embedded twin and α (x ) = −α U ˜ is given by: in the untwinned region. Accordingly, in view of Eqs. (17)–(18), the discontinuous operator B
δik δ jl + α Z Ti j Tkl B¯ klA in the embedded twin, ˜ i jA = B ˜U B = δik δ jl − α U Ti j Tkl B¯ klA in the untwinned region, i jA
˜ Zi jA = B
(31)
¯ and unit normal to the twinning plane n, and the unknown parameters m, in terms of its known continuous counterpart B α Z and α U , which are determined as described below. Integrals over the element domain e , such as the one in Eq. (28), are approximated via Gaussian quadrature as follows:
e
f ( x ) d ≈
Nq
wq J ( ξ q ) f ( ξ q ),
(32)
q=1
where Nq is the number of integration points, ξ q and wq are the coordinates of the q-th integration point in e and its corresponding weight coefficient, respectively, and J is the Jacobian determinant of the isoparametric mapping. After the twin is embedded within an element, the integrand f(ξ q ) assumes two different values at each integration point, fZ (ξ q ) and fU (ξ q ), representing the embedded twin and the surrounding untwinned material, respectively. The following modified quadrature rule, introduced by Belytschko et al. (1988) and Fish and Belytschko (1988) , is adopted in this case to evaluate the total integral:
e
f ( x ) d ≈
Nq
wq J ( ξ q )
q=1
b Z f (ξ q ) + h
1−
b h
f U (ξ q ) ,
(33)
√ where b is the width of the embedded twin, and h = A is the characteristic size of the element e , whose area is denoted by A. Based on the above modified quadrature rule, the orthogonality condition (27) can be expressed as:
e
¯ −B ˜ d ≈ B
Nq q=1
wq J (ξ q ) B¯ (ξ q ) −
b Z ˜ (ξ q ) + B h
1−
b ˜ U (ξ q ) B h
= 0.
(34)
Substituting (31) into (34), we obtain the following expression for the parameter α Z in terms of α U :
α = Z
h −1 b
αU.
(35)
8
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Similarly, based on Eq. (28), the element internal force vector is written as follows:
fAint
e
=
˜ i jA σi j d ≈ B
e
Nq
wq J ( ξ q )
q=1
b ˜ Zi jA (ξ q ) σiZj (ξ q ) + B h
1−
b U ˜U B i jA (ξ q ) σi j (ξ q ) , h
(36)
where σ Z and σ U denote the Cauchy stress in the twin and the untwinned region, respectively. Noting that the same slipbased mechanisms govern the evolution of plastic deformation in these two regions, the stresses σ Z and σ U are computed using the same single-crystal plasticity model described in Section 2.3. It is important to note, however, that the slip systems α {Sˆ 0 } and {Sα0 }, on which dislocation motion takes place in the twinned and untwinned regions, respectively, have different orientations, and the plastic deformation rates, LZ,p and LU,p , in these two regions are different as well; cf. Eqs. (8) and (11). Equilibrium at the interface between the embedded twin and the surrounding untwinned material requires σ Z and σ U to satisfy the following traction-continuity condition:
n ·
Z Z,e Z,p Z,p σ F , F , L − σ U FU,e , FU,p , LU,p d = 0
(37)
where d is an element of area on the twin boundary in the deformed configuration. The above nonlinear equation is solved for m and α U using the Newton iterative scheme detailed in Appendix A. 2.2.5. Time integration The time interval of interest, T , is split into increments [t n , t n+1 ] of size t := t n+1 − t n . Quantities evaluated at the beginning and end of a typical time increment are denoted, respectively, by (•)n and (• )n+1 . Each time step, the nodal accelerations d¨ n+1 are obtained by solving the dynamic finite element equations (21) explicitly: −1
ˆ d¨ n+1 = M
int f ext n+1 − f n+1 ,
(38)
ˆ for improved computational efficiency, and then the replacing the consistent mass matrix M by its lumped counterpart M central difference scheme is used to update the nodal solution variables as follows:
1 d¨ n t 2 , 2
(39)
1 d¨ n + d¨ n+1 t. 2
(40)
dn+1 = dn + d˙ n t + d˙ n+1 = d˙ n +
The global internal force vector f int n+1 in Eq. (38) is computed via standard finite element assembly procedures (Hughes, 20 0 0) from its elemental counterpart given by Eq. (36) in terms of the updated stresses σ Zn+1 inside the twin and σ Un+1 in the untwinned material. These stresses are in turn computed using the single-crystal elasto-viscoplasticity model detailed in several published works (Anand, 2004; Kalidindi et al., 1992; Kalidindi, 1998) and outlined in Section 2.3. This constitutive model requires, as input, a deformation gradient that consists of elastic and plastic (slip) parts only. The appropriate deformation gradient tensors are obtained as follows: Recasting Eq. (3) as F˙ = LF, and using the mid-point exponential map integrator, we obtain
Fn+1 = exp[ L]Fn ≈ [1 + L]Fn ,
(41)
where
L := tLn+ 12 = t
∂ vn+ 12 ∂ u ∂ (un+1 − un ) = = . ∂ xn+ 12 ∂ xn+ 12 ∂ xn+ 12
(42)
As noted above, in the present finite element formulation, the stress must be evaluated using the constitutive model as a function of the assumed velocity gradient L. Therefore, based on (41), we write
Fn+1 ≈ [1 + L]Fn
⇒
Fn+1 =
FZn+1 ≈ 1 + LZ FZn FU n+1
≈ 1 + L
U
FU n
in the embedded twin, in the untwinned region,
(43)
where, in light of Eq. (30),
Z Li j = δik δ jl + α Z Ti j Tkl Lkl in the embedded twin, L i j = LUi j = δik δ jl − α U Ti j Tkl Lkl in the untwinned region.
(44)
The stress in the untwinned material, σ U , is computed using the single-crystal plasticity model as a function of FU , n+1 n+1 which can be decomposed into elastic and plastic (slip) parts according to Eq. (7). On the other hand, the stress within the twin, σ Zn+1 , is computed using the same material model as a function of
Fˆ Zn+1 := FZn+1 FZ,tw
−1
,
which can be decomposed multiplicatively into elastic and plastic (slip) parts according to Eq. (9).
(45)
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
9
2.3. Constitutive modeling Based on the kinematic relationships described in Section 2.1 and the assumption that thermally activated slip is the dominant mechanism for plastic deformation in the twinned and untwinned regions alike, the stress in these two regions is computed using the same single-crystal elasto-viscoplasticity model. This constitutive model is outlined below, and a more detailed description thereof can be found in the literature (Anand, 2004; Kalidindi et al., 1992; Kalidindi, 1998). Although the model includes temperature as a state variable, all simulations presented herein are conducted assuming isothermal conditions. For convenience, we employ notation that does not differentiate between quantities associated with the twinned or untwinned regions in the following brief description of the material model (Section 2.3.1) and its algorithmic implementation (Section 2.3.2). These differences are highlighted in Section 2.3.3, and the material Jacobian is discussed in Section 2.3.4.
2.3.1. Single-crystal elasto-viscoplastic model The theoretical foundations of single-crystal elasto-viscoplasticity can be traced back to the work of Asaro and Needleman (1985) Asaro and Rice (1977), Asaro (1983), Hill and Rice (1972), and Rice (1971). Let Cc represent the elasticity tensor in the crystal basis that linearly depends on the temperature θ as follows:
Cci jkl (θ ) = C0i jkl + Mi jkl θ ,
(46)
where C0i jkl is the elasticity tensor at θ = 0 K, and Mi jkl are constant coefficients. Let Q represent the orthogonal tensor that rotates the crystal basis to the global basis. Then, the fourth-order elasticity tensor in the global basis is given by
Ci jkl = Qim Q jn Qkp Qlq Ccmnpq .
(47)
The free energy is expressed as follows:
1 2
ψ ( Ee ) = Ee : C ( θ ) : Ee ,
(48)
where Ee is the elastic Green–Lagrange strain tensor defined in the relaxed (second) intermediate configuration as:
Ee =
1 e 1 ( C − 1 ) = Fe T Fe − 1 . 2 2
(49)
The second Piola–Kirchhoff (2nd P-K) stress Se in the (second) intermediate configuration, which is work-conjugate to the elastic strain tensor Ee , is calculated as:
Se = C : Ee .
(50)
Recall that the plastic deformation rate Lp in Eq. (4) is based on the Schmid tensor and the plastic slip rate. The flow rule for the slip rate γ˙ α on the α -th slip system is expressed as (Bronkhorst et al., 2007; Busso, 1990):
F γ˙ α = γ˙ 0 exp − 0 1 − kθ
|τ α | − sαρ μ/μ0 sαl μ/μ0
p q
sign(τ α ),
(51)
where k is the Boltzmann constant, γ˙ 0 is a constant reference rate of shearing, sα ρ is the deformation resistance due to evolving dislocation density, sα is the constant intrinsic lattice resistance, the Macaulay brackets • designate the ramp l function, i.e. x := max(x, 0 ), ∀x ∈ R, the shear modulus ratio μ(θ )/μ0 encapsulates the temperature sensitivity, and τ α is the resolved stress on the α -th slip system, defined as:
τ α = (Ce Se ) : Sα0 .
(52)
The evolution of the deformation resistance sα ρ on slip system α is described as:
s˙ αρ =
hαβ |γ˙ β |,
hαβ = r + (1 − r )δαβ hβ ,
β
β
β
ss − sρ hβ = h0 β β ss − s0
(53)
where r is the latent-hardening parameter (Peirce et al., 1982), hβ is the self-hardening rate including both the dislocation β generation and annihilation effects (Acharya and Beaudoin, 20 0 0), and ss is the saturation stress parameter that depends on the temperature and shear rate as follows (Kocks, 1976): β
β
ss = ss0
|γ˙ β | γ˙ 0
kAθ .
(54)
10
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
2.3.2. Algorithmic implementation We use the compact notation t = tn and τ = tn+1 to represent, respectively, the beginning and the end of a typical time interval [tn , tn+1 ]. Letting γ α = γ˙ α (τ ) t represent the increment of the plastic shearing on the α -th slip system, and using Eq. (4) we write:
F (τ ) = exp t p
α
γ˙ α (τ )Sα0 Fp (t ) = exp
Therefore,
Fp (τ )−1 = Fp (t )−1 exp −
α
α
γ α Sα0 Fp (t ) ≈
1+
α
γ α Sα0 Fp (t ).
γ α Sα0 ≈ Fp (t )−1 1 − γ α Sα0 .
(55)
(56)
α
Based on the multiplicative decomposition, F(τ ) = Fe (τ )Fp (τ ), we have:
F ( τ ) = F ( τ )F ( τ ) e
p
−1
≈ F(τ )F (t ) p
−1
1−
α
γ α Sα0 .
(57)
Assume that the deformation increment F = F(τ )F(t )−1 in the time step is purely elastic, we define the trial elastic deformation at the end of the time step as follows:
Fe trial (τ ) = FFe (t ) = FF(t )Fp (t )−1 = F(τ )Fp (t )−1 .
(58)
As a result, the relationship between Fe (τ ) and Fe trial (τ ) can be expressed as:
Fe (τ ) ≈ Fe trial (τ ) 1 −
α
γ α Sα0 .
(59)
Based on Eq. (59), the right Cauchy–Green tensor is expressed as:
C (τ ) = F (τ ) F (τ ) ≈ e
e
T e
1−
=
1−
α
α
γ α (Sα0 )T Fe trial (τ )T Fe trial (τ ) 1 −
α
γ α Sα 0
γ α (Sα0 )T Ce trial (τ ) 1 − γ α Sα0
≈ Ce trial (τ ) −
α
α
γ α Ce trial (τ )Sα0 + Sα0 T Ce trial (τ ) ,
(60)
where the second-order term O (( γ α )2 ) is neglected. Therefore, the elastic strain is written as:
Ee (τ ) ≈ Ee trial (τ ) −
α
γ α sym Ce trial (τ )Sα0 ,
(61)
and the 2nd P-K stress Se (τ ) at the end of the time step is given by:
Se (τ ) = C : Ee (τ ) ≈ C : Ee trial (τ ) − = Se trial (τ ) −
α
α
γ α sym Ce trial (τ )Sα0
γ α C α ( τ ) ,
(62)
where Cα (τ ) = C : sym[Ce trial (τ )Sα ]. 0 α In a single crystal, the prescribed slip systems in the reference configuration {t α 0 , s0 }, {α = 1, 2, . . . , N }, are timeindependent. In the time interval [t, τ ], the deformation gradients at the beginning and the end of the time step, F(t) and F(τ ), are known. In a finite element simulation, the role of the constitutive update subroutine is to compute the variables {Fp (τ ), sαρ (τ ), Se (τ )} at the end of the time step, starting from their known values {Fp (t ), sαρ (t ), Se (t )} at the beginning of the time step. This update procedure is governed by a set of implicit equations, summarized as follows, solved using a two-level iterative procedure detailed by Anand (2004) and Kalidindi et al. (1992):
Se (τ ) = Se trial (τ ) − α γ α Cα (τ ), sαρ (τ ) = sαρ (t ) + β hαβ (τ )| γ β |, γ α = γ˙ α (τ ) t.
(63)
2.3.3. Differences between twinned and untwinned regions For conciseness, a notation that does not differentiate between quantities associated with the twinned and untwinned regions has been employed in the foregoing overview of the constitutive equations and stress update procedure. For the sake of clarity, several important points are emphasized here:
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
11
1. The tensors Ee and Se are defined on the intermediate configuration in the untwinned region, and on the second intermediate configuration in the twinned region; cf. Figs. 1 and 2. α 2. When the resolved shear stress τ α is computed according to Eq. (52), the Schmid tensors {Sα } and {Sˆ 0 } are used in 0 the untwinned and twinned regions, respectively. 3. The elasticity tensors CU and CZ corresponding to the untwinned and twinned regions, respectively, are related as follows (in the global basis):
CZi jkl = Rim R jn Rkp Rlq CU mnpq ,
(64)
where R is the orthogonal matrix, defined in Eq. (13), which relates the lattice orientations in the twinned and untwinned regions of the crystal. 4. In the stress update procedure, FU = FU,e FU,p is provided as input to compute the stress σ U in the untwinned region, while Fˆ Z := FZ (FZ,tw )−1 = FZ,e FZ,p is provided as input to compute the stress σ Z inside the twin. 2.3.4. Material Jacobian In order to enforce the traction continuity condition (37), the Newton scheme described in Appendix A requires the calculation of the material Jacobian W (τ ) based on the constitutive equations summarized in the present section,
W (τ ) =
∂ σ (τ ) , ∂ F (τ )
(65)
where σ is the Cauchy stress tensor obtained by pushing forward the 2nd P-K stress Se as follows:
σ=
1 T Fe Se Fe . detFe
(66)
The derivation of the material Jacobian W is presented in Appendix B. Note that the Jacobian derived by Balasubramanian (1998) and Dai (1997), is defined as follows:
¯ := W
∂σ ∂σ ≈ , ∂E ∂U
based on the approximation
E = lnU ≈ U − 1, where U is the right stretch tensor obtained via polar decomposition of the deformation gradient, F = RU. In Appendix B, we derive the expression for the Jacobian defined in Eq. (65), which is required in the embedded-discontinuity finite element formulation. 2.4. Twin nucleation, propagation and growth models The twin initiation criteria adopted in this paper is based on a probabilistic twin nucleation model detailed in (Beyerlein and Tomé, 2010; Niezgoda et al., 2014). One of the core assumptions of this probability-based model is that twin nucleation occurs when some number of grain boundary dislocations undergo stress-driven transformations that then coalesce into a single stable nucleus. The cumulative total number of such transformations N(τ ) on the stress interval [0, τ ] is treated as a random variable and described by a Poisson process. The probability of observing N = k transformations on the stress interval [0, τ ] is expressed as:
P (N (τ ) = k ) =
(τ )k k!
e−(τ ) .
(67)
In Eq. (67), (τ ) represents the cumulative rate of the grain boundary dislocation transformations, defined as:
(τ1 , ) =
τ1 + τ1
λ ( τ ) dτ ,
(68)
where is the stress interval, and λ(τ ) represents the probability of the grain boundary dislocation transformations. In order to generate a stable twin nucleus, it is assumed that at least k∗ transformation events are required and this reorganization is a stress driven process. Let the random variable S represent the critical nucleation strength, or the stress required to transform an appropriate number of grain boundary dislocations, then,
P (τ > S ) = P (N (τ ) ≥ k ) = 1 − ∗
∗ k −1
P ( N ( τ ) = k ).
(69)
k=0
Based on the assumption that k∗ = 1 and all dislocation events generate a stable twin nucleus, we have
P (τ > S ) = P (N (τ ) ≥ 1 ) = 1 − e−(τ ) .
(70)
12
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
In order to decide the specific form of the cumulative rate of the grain boundary dislocation transformations (τ ), the corresponding probability λ(τ ) is assumed to follow the Weibull distribution, so that λ(τ ) → 0 as τ → 0, preventing the case of non-zero transformation rate at zero stress. As a result,
(τ ) =
τ β τc
,
(71)
where τ c is the Weibull scale parameter and can be interpreted as a characteristic stress for nucleation, and β is the Weibull shape parameter. Both τ c and β are treated as material constants and do not depend on temperature or loading conditions. In the context of single-crystal elasto-viscoplasticity presented in Section 2.3, τ is the resolved shear stress associated with α ˆ 0 under consideration. In the finite element simulation of twin propagation using the embedded weak the twin system S discontinuities introduced in Section 2.2, once the twin nucleation probability associated with an element, calculated according to Eq. (70), exceeds the threshold value P0 , that is,
P (τ > S ) = 1 − e−(τ ) > P0 ,
(72)
we assume that the twin has nucleated inside this element. In the finite element simulation, all the state variables, including the twin nucleation probability, are evaluated at the Gauss quadrature points. In the current implementation, the twin nucleation probability within an element is taken as the highest value calculated at the Gauss points inside the element. An alternative approach is to evaluate Eq. (72) at the element centroid. After twin nucleation is detected at a grain boundary, the generated twin starts to grow, and eventually propagates into neighboring elements under subsequent deformation. Let ltw represent the embedded twin length inside individual elements, and recall that b represents the embedded twin thickness. We adopt the following expressions for the twin propagation rate l˙tw and the twin thickness growth rate b˙ tw :
−Gtw l˙tw = vl exp 1− kθ
τ p q sl
,
−Gtw b˙ tw = vb exp 1− kθ
τ p q sl
,
(73)
where vl is the reference propagation speed, and vb is the reference thickness growth speed. During a typical time interval [tn , tn+1 ], the embedded twin length and the twin thickness at the end of the time step are calculated as: tw ˙tw lntw +1 = ln + ln+1 t,
bn+1 = bn + b˙ tw n+1 t.
(74)
Once the twin length ltw inside an element is larger than the length of the line segment le bounded by the intersection points between the twin plane and the element boundary, that is,
l tw ≥ l e ,
(75)
the twin is allowed to propagate into the neighboring element. We summarize the computational framework, including the twin initiation and propagation scheme, in Algorithm 1 , where we use the following terminology and notation. The grain boundary set, denoted by Egb , is the set of elements located on the grain boundaries between different neighboring crystals. The embedded twin set, denoted by Etw , is the set of e , is the set of neighboring elements elements containing the embedded twin. The neighboring element set, denoted by Enb e e that share a common edge with an individual element . The sets of Egb and Enb are known at the start of the finite element simulation from the initial arrangement of multi-crystals and the spatial discretization, and the set of Etw is empty before twin nucleation is detected. 3. Numerical results In this section, we demonstrate the capability of the proposed computational framework through several numerical examples in two dimensions (2D). Specifically, we focus on simulating the initiation and propagation of the {101¯ 2} tensile twin system in Titanium. In the first numerical example, we construct a simple tri-crystal structure, and assume that a fully developed {101¯ 2} tensile twin is embedded in the center of the middle crystal (Grain C) from the beginning of the simulation. We aim to verify the embedded weak discontinuity formulation by comparing the corresponding numerical results with the ones obtained from the conventional finite element simulation that explicitly resolves the embedded tensile twin. In the second numerical example, we use the same tri-crystal structure. However, instead of assuming that the location of the twin is known a priori, we rely on the probability-based twin nucleation criterion presented in Section 2.4 to detect the initiation of the twin. Moreover, instead of assuming that the twin is fully developed from the beginning, we allow the twin to evolve (propagate and grow) according to the growth rates expressed in Eq. (73). In the third numerical example, we simulate the initiation and propagation of the {101¯ 2} tensile twin in a more general poly-crystalline structure. The material parameters used for HCP titanium are given in Table 1. A constant temperature of 295 K is assumed everywhere in the domain in all the simulations. The simulations for all the aforementioned problems are conducted in two-dimensions. In reality, however, structures of crystal lattices, and their slip and twin systems, are three-dimensional. Therefore, spatial rotations are necessary to simplify the {101¯ 2} twin propagation that is intrinsically three-dimensional, to a two-dimensional setup. Let Tc represent a secondorder tensor defined in the orthogonal basis {eci } associated with the crystal lattice, and T represent the counterpart defined
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
13
Algorithm 1: Main steps performed within a single time increment [tn , tn+1 ].
1
2 3
• Loop over elements in the discretized domain h : if e ∈ Egb (current element is located at the grain boundary) then Calculate the twin nucleation probability Pi according to Eq. (70) at i-th Gauss point in e ; if max {Pi } > P0 (P0 is the threshold value) then Add e to the embedded twin set Etw ; Compute the twin plane normal in the current configuration via push-forward of its (known) counterpart in the undeformed configuration: n =
4 5
6
7 8 9 10 11
12
F −T n n+1 0
;
F −T n n+1 0
Store the coordinate of the element centroid as x0 ; Obtain the twin propagation trajectory {x|(x − x0 ) · n = 0}; if e ∈ Etw (current element has embedded twin) then , α Z and α U by solving the traction continuity condition (37) using the Newton scheme detailed in Calculate m Appendix A; Z and σ U , in the twinned and untwinned regions, respectively; Update stresses σ Compute the contribution of e to the internal force fAint e , using Eq. (36); Calculate the twin propagation rate l˙itw and the thickness growth rate b˙ tw according to Eq. (73) at i-th Gauss point; i Update the twin length litw and the twin thickness btw inside the element according to Eq. (74) at i-th Gauss point; i Calculate the length of the line segment l e bounded by the intersection points between the twin propagation trajectory and the boundary of the element, {x|(x − x0 ) · n = 0} ∩ ∂ e ; if max {litw } ≥ l e then
• Loop over the neighboring elements enb ∈ Enb of the current element e :
if enb ∈ / Etw and x|x ∈ enb ∩ {x|(x − x0 ) · n = 0} = ∅ (the neighbor element enb intersects with the twin
13
14 15
propagation trajectory) then Add enb to the embedded twin set Etw ; else U in the untwinned material at all Gauss points in e ; Update stress σ Compute the contribution of e to the internal force fAint e , using Eq. (36) with b = 0;
Table 1 Material parameters used for HCP titanium. Parameter
Value
Unit
Parameter
Value
Unit
ρ θ
4507 298.0 177.59 191.34 87.57 45.01 52.23 −52.0 −36.7 13.3 −32.8 −19.7 1.0 × 107 250.0
kg/m3 K GPa GPa GPa GPa GPa GPa/K GPa/K GPa/K GPa/K GPa/K s−1 MPa
F0 p q r h0 s0 A ss0
3.0 × 10−19 0.6 1.4 1.4 225.0 30.0 1.0 × 10−19 400.0 11.1 200.0 0.85 3190.0 10.0 2.76 × 10−19
J — — — MPa MPa J MPa — MPa — m/s m/s J
C011 , C022 C033 C012 , C013 , C023 C044 C055 , C066 M11 , M22 , M33 M12 , M13 , M23 , M44 M55 , M66 γ˙ 0 sl
β τc
P0 vl vb Gtw
in the global Cartesian coordinate basis {ei }. The two tensors T and Tc are related through a rotation matrix Q in the following form:
T = Q Tc QT . The rotation matrix Q rotates the crystal basis c
e = Qe .
(76)
{eci }
to the global basis {ei } according to the following relationship:
(77)
The rotation represented by the matrix Q can be decomposed into a (non-unique) sequence of three simple rotations using Euler angles {φ , θ , ω}. The three simple counter-clockwise rotations adopted in this paper are enumerated as follows (Kalidindi et al., 1992):
14
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Fig. 4. The tri-crystal problem, where the fully developed {101¯ 2} tensile twin is embedded in the center of Grain C at the beginning of the simulation: (a) the boundary conditions, (b) the velocity profile applied at the top surface (corresponding to a shear strain rate of 7.5 × 104 s−1 ).
Fig. 5. Rotation of the crystal lattice so that, initially (at t = 0), the plane spanned by the normal and shear directions (n and m) of the {101¯ 2} tensile twin system coincides with the e1 –e2 plane (i.e. the x–y plane) in the global coordinate system.
(1) Given an orthogonal basis {x, y, z} aligned with {ec1 , ec2 , ec3 } associated with the crystal lattice, rotate {x, y, z} by the first Euler angle φ about the z-axis to get {x(1) , y(1) , z}. This is represented by:
Q1 =
cos φ −sin φ 0
sin φ cos φ 0
0 0 . 1
(2) Rotate {x(1) , y(1) , z} by the second Euler angle θ about the x(1) -axis to get {x(1) , y(2) , z(1) } so that z(1) is aligned with the out-of-plane direction e3 . This is represented by:
Q2 =
1 0 0
0 cos θ −sin θ
0 sin θ . cos θ
(3) Rotate {x(1) , y(2) , z(1) } by the third Euler angle ω about the z(1) -axis to get {x(2) , y(3) , z(1) } so that it is aligned with {e1 , e2 , e3 } of the global basis. This is represented by:
Q3 =
cos ω −sin ω 0
sin ω cos ω 0
0 0 . 1
Combining the three simple rotations described above, the rotation matrix Q is expressed as:
Q = Q3 Q2 Q1 =
cos φ cos ω − sin φ cos θ sin ω −cos φ sin ω − sin φ cos θ cos ω sin φ sin θ
sin φ cos ω + cos φ cos θ sin ω −sin φ sin ω + cos φ cos θ cos ω −cos φ sin θ
sin θ sin ω sin θ cos ω . cos θ
(78)
3.1. Tri-crystal with a pre-existing twin We construct the tri-crystal problem shown in Fig. 4, where three grains with different orientations undergo a simple shear in the x-direction. The tri-crystal has a total width of 200 μm and a height of 50 μm. Grains A and B have a width of 50 μm each. Fig. 4(a) shows the boundary conditions used in the simulation, where the horizontal velocity profile shown in Fig. 4(b) is applied to the top surface of the tri-crystal sample, corresponding to a shear strain rate of 7.5 × 104 s−1 . Based
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
15
Fig. 6. With the conventional FEM, (a) the {101¯ 2} tensile twin is explicitly resolved using a refined finite element mesh. With the embedded weak discontinuity approach, (b) a coarser mesh is used in the parent grains/untwinned region, and (c) the twin is represented by weak discontinuities embedded within a subset of the elements (highlighted in blue). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 7. Contours of the shear stress σ 12 (unit: MPa) at very early stages in the simulation, showing as the stress wave (a) begins propagating from the top surface, (b) reaches the twinned region, (c) reaches the bottom surface, and then (d) bounces back and starts propagating upward.
on the rotation matrix defined in Eq. (78), the Euler angles {φ = 210◦ , θ = 90◦ , ω = 42.472◦ } are used to rotate the crystal lattice associated with Grain C (the middle grain), as shown in Fig. 5, so that the plane spanned by the normal and shear directions {n, m} of the {101¯ 2} tensile twin system in Grain C coincides with the x–y plane in the global coordinate system. Thus, before the tri-crystal deforms (at t = 0), the normal to the twin plane points in the y-direction, i.e. n = [0, 1, 0]T , and the shear direction points in the x-direction, i.e. m = [1, 0, 0]T , in the global coordinate system. Grains A and B are both assigned the Euler angles {φ = 210◦ , θ = 90◦ , ω = 27.472◦ }, so that two 15◦ tilt grain boundaries are formed (see Fig. 4). With the aid of these crystal rotations, the problem is cast in a two-dimensional (2D) setting. We assume that the {101¯ 2} tensile twin has fully developed inside Grain C from the beginning, as shown in Fig. 4, and the location of the twin is known a priori (in the middle of Grain C). Based on the transformation matrix between the lattice orientations in the untwinned and the twinned regions, as shown in Eq. (13), the Euler angles of the crystal lattice inside the embedded twin is {φ = 30◦ , θ = 90◦ , ω = −42.472◦ }. In order to verify the proposed computational framework, we use two different approaches in the numerical simulations. In the first approach, we use a refined mesh to explicitly resolve the tensile twin, as shown in Fig. 6(a). The numerical results obtained using this approach are treated as the reference solution. In the second approach, we adopt a coarser mesh and use the proposed computational framework to treat the tensile twin as embedded weak discontinuities inside Grain C, as shown in Figs. 6(b) and (c). Two cases are considered. In the first case, the thickness b of the twinned region is 4.5 μm, and thus the ratio between the twin thickness and the element characteristic length is close to unity (b/h ≈ 1.0). Fig. 7 shows the propagation of the
16
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Fig. 8. Distributions of the third Euler angles ω (unit: degree) in Grains A and B, and in the twinned and untwinned regions of Grain C. Results obtained using the embedded weak discontinuity appraoch with b/h ≈ 1.0 (b and c) are in very good agreement with those obtained using the conventional FEM (a).
Fig. 9. Comparison of the conventional FEM and the embedded weak discontinuity simulations: (a) evolution of the third Euler angle ω in the twinned and untwinned regions, and (b) the load–displacement curves.
shear stress σ 12 in very early stages of the simulation conducted using the proposed computational framework. As expected, once the velocity profile shown in Fig. 4(b) is applied to the top surface, the stress wave propagates downward until it reaches the bottom surface, and then bounces back. Fig. 8 shows the distributions of the third Euler angle ω in the different grains at the end (t = 10−6 s) of the simulations conducted using the conventional finite element method and the embedded weak discontinuity approach. The lattice rotations in the twinned and untwinned regions obtained using the two numerical approaches are in very good agreement, which is confirmed by the evolution history of the third Euler angle ω reported in Fig. 9(a). Moreover, the load–displacement plots generated using the embedded weak discontinuity and conventional FEM results, are indistinguishable as can be seen from Fig. 9(b). Fig. 10 further examines other state variables, including the cumulative slip and the stress components inside the twinned and untwinned regions, obtained using the two numerical approaches. Again, the numerical results obtained using the embedded weak discontinuity approach, shown in Fig. 10(b), (c), (e), (f), (h), (i), (k) and (l), are in good agreement with the corresponding results obtained using the conventional FEM with a refined mesh that explicitly resolves the twin, shown in Fig. 10(a), (d), (g) and (j). Note that among the three stress components σ 11 , σ 22 and σ 12 , only the first one (σ 11 ) is different in the twinned region (Fig 10(f)) compared to the untwinned region (Fig. 10(e)). The other two stress components (σ 22 and σ 12 ) have the same magnitudes in the twinned and untwinned regions; see Fig. 10(h), (i), (k) and (l). This is due to the enforcement of the traction continuity condition (37), along the (horizontal) interface between the twinned and untwinned regions. In the second case, the thickness b of the fully developed twin is 1.5 μm, and the ratio between the twin thickness and the element characteristic length is b/h ≈ 1/3. Figs. 11 and 12 compare the results obtained using the conventional FEM and
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
17
Fig. 10. Distributions of the cumulative slip on all slip systems (a–c) and the in-plane stress components (d–l, unit: MPa) obtained using the conventional FEM and the embedded weak discontinuity approach with b/h ≈ 1.0. The results obtained using the proposed framework (b, c, e, f, h, i, k, l) are in very good agreement with the corresponding results obtained using the conventional FEM (a, d, g, j).
Fig. 11. Distributions of the 3rd Euler angles ω (unit: degree) in Grains A, B, and the twinned and untwinned regions of Grain C. The numerical results obtained using the proposed framework with b/h ≈ 1/3 (c–f) are in good agreement with the corresponding results obtained using conventional FEM (a–b). Note that the band of elements highlighted in blue in (e–f) is about three times wider than the twinned region embedded within these elements. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
18
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Fig. 12. Distributions of the cumulative slip on all slip systems (a–c) and the in-plane stress components (d–l, unit: MPa) obtained using the conventional FEM and the embedded weak discontinuity approach with b/h ≈ 1/3. The numerical results obtained using the proposed framework (b, c, e, f, h, i, k, l) are in good agreement with the corresponding results obtained using conventional FEM (a, d, g, j). Note that the band of elements highlighted in blue in (e–f) is about three times wider than the twinned region embedded within these elements. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 13. Twin nucleation probability distribution in the elements adjacent to the grain boundaries (a) at t = 0.60 × 10−6 s, and (b) at t = 0.646 × 10−6 s.
the embedded weak discontinuity approach. Note that in these two figures, distributions of state variables in the twinned region are plotted over the whole area of the elements in which the twin is embedded, even though the twinned region only occupies one third of the area of those elements. As in the first case, we observe good agreement between the results of the embedded weak discontinuity approach and the conventional FEM simulation for the naive tri-crystal problem, in which the fully developed twin is embedded into the mesh at the beginning of the simulation. The above numerical simulations are used as a simple test to verify the implementation of the proposed computational framework.
3.2. Twin initiation and propagation in tri-crystals In this numerical example, we use the same geometry and boundary conditions of the previous tri-crystal problem, shown in Fig. 4. However, in this case, we adopt the probability-based twin nucleation model, introduced in Section 2.4, to detect the location of twin initiation at the grain boundaries. Moreover, instead of assuming that the twin is fully developed from the beginning, we allow the twin to evolve according to the evolution Eq. (73). That is, after twin initiation is detected at one grain boundary, the twin thickness and length both grow under subsequent deformation (b and ltw are functions of time t), and the twin propagates into neighboring elements once inequality (75) is satisfied. Fig. 13 shows the distribution of the twin nucleation probability at the two grain boundaries. Once initiation of the {101¯ 2} tensile twin is detected at one of the grain boundaries (at t = 0.323 × 10−6 s), the twin starts to grow in thickness and propagate along the twinning plane into the neighboring grain. Fig. 14 shows the twin propagation and the growth of the twin thickness b(t) during the deformation process. After the twin is initiated at the left grain boundary, it starts to
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
19
Fig. 14. Growth and propagation of the {101¯ 2} tensile twin during the deformation process in the tri-crystal problem. Note that the characteristic element size in this problem is h ≈ 5.7 μm.
Fig. 15. Stress components σ 11 and σ 12 (unit: MPa) in the untwinned (a, c) and twinned (b, d) regions at t = 0.660 × 10−6 s.
thicken while at the same time propagating to the right. This propagation process continues until the twin reaches the right grain boundary. Fig. 15 shows the distributions of the stress components σ 11 and σ 12 inside the untwinned and twinned regions at t = 0.660 × 10−6 s, approximately when the twin reaches the mid-point between the left and right grain boundaries. During this process, the traction continuity condition n · [[σ ]] = 0 is enforced along the horizontal interface between the twin and
20
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Fig. 16. Evolution of the stress components σ 11 and σ 12 (unit: MPa) inside the untwinned (left) and twinned (right) regions after the twin is fully developed: The σ 11 component exhibits stress relaxation in the embedded twin (b, f, j) but remains almost constant in the untwinned region (a, e, i); meanwhile, the σ 12 component remains equilibrated ([[σ12 ]] = 0) between the twinnned and untwinned material, due to proper enforcement of the traction continuity condition n · [[σ ]] = 0.
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
21
Fig. 17. Poly-crystal problem: (a) the geometry and boundary conditions, (b) the velocity profile applied at the top surface, corresponding to a strain rate of 15.0 × 105 s−1 , and (c,d) Finite element mesh and value of the third Euler angle ω in each grain at t = 0.
the surrounding material. This can be confirmed, for example, by comparing Fig. 15(c) and (d), and observing that σ 12 has the same value in the twinned region and the surrounding untwinned material (i.e. the jump [[σ12 ]] = 0). Conversely, since the interface is horizontal (i.e. n points in the x2 -direction), traction continuity does not impose any constraint on the σ 11 component, and thus the jump [[σ11 ]] = 0, as can be seen by comparing Fig. 15(a) and (b). Fig. 16 shows the distributions of the same two stress components inside the untwinned and twinned regions at three different points in time, and shows that the same observations hold true after the conclusion of the twin propagation stage. It is also observed by comparing Fig. 15(a) to (b), and Fig. 16(a) to (b), that initially, σ 11 is larger inside the twin than in the surrounding untwinned material. However, stress relaxation takes place very shortly thereafter inside the twinned region, as can be seen by comparing Fig. 16(b), (f), and (j). This stress relaxation is due to dislocation slip-based plastic deformation, which is concentrated inside the twinned region as can be seen by comparing Fig. 10(b) to (c), and Fig. 12(b) to (c). The short delay associated with this effect can be attributed to the rate-dependent nature of the plastic response, represented by the single crystal constitutive model described in Section 2.3.1.
22
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Fig. 18. Twin nucleation probability distributions in the elements adjacent to the grain boundaries at different time steps. Two twin nucleation events are detected. Results obtained using a mesh size h = 2 μm (a–c) are compared to those obtained using a coarser mesh with h = 4 μm (d–f).
3.3. Twin initiation and propagation in poly-crystal In this example, the numerical simulation is performed on a 100 × 100 μm2 poly-crystalline structure that undergoes uniaxial tension, as shown in Fig 17. Figs. 17(a) and (b), respectively, show the boundary conditions and the velocity profile applied at the top surface, corresponding to a strain rate of 15 × 105 s−1 . A constant temperature of 295 K is assumed everywhere in the domain. Two different meshes, with the average mesh size h = 4 μm and h = 2 μm, are used to discretize the poly-crystalline domain, as shown in Fig. 17(c) and (d). Each grain is assigned a different third Euler angle ω in order to form tilt grain boundaries. Fig. 18 compares the distributions of twin nucleation probability at the grain boundaries obtained from simulations using different mesh sizes. Fig. 18(a), (b) and (c) show the results obtained at different time steps from the simulation with the mesh size h = 2 μm, and Figs. 18(d), (e) and (f) show the results obtained at different time steps from the simulation with the mesh size h = 4 μm. Fig. 18 shows that two twin nucleation events are detected at the grain boundaries. Due to the influence of the mesh size h on the magnitude of the stress calculated at the grain boundary, the second twin nucleation event is detected at slightly different time values, depending on the mesh size used in the simulation. However, by comparing Fig. 18(b) with (e) and Fig. 18(c) with (f), we can see that the first and second twin nucleation sites are mesh-independent. After twin initiation is detected, the two sets of twins start to propagate and grow in thickness. Fig. 19 shows the propagation paths of the tensile twins initiated at the grain boundaries. In this figure, the dark gray area represents elements adjacent to the grain boundaries, and the red area represents the elements containing the embedded twin. It is worth emphasizing that the size of the elements highlighted in red does not represent the thickness of the embedded twin. Once initiated from the grain boundaries, the two tensile twins start to propagate into the middle grain, and the propagation continues until they reach another grain boundary. We also compare the twin propagation paths obtained from the simulations using different mesh sizes. Fig. 19(a), (b) and (c) are obtained from the simulation using the mesh size h = 2 μm, and
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
23
Fig. 19. Propagation paths of the two twins initiated at the grain boundaries. The red areas consist of elements containing an embedded twin, and the dark gray areas comprise elements that are adjacent to the grain boundaries. Results obtained using a mesh size h = 2 μm (a–c) are compared to those obtained using the coarser mesh with h = 4 μm (d–f). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 19(d), (e) and (f) are obtained from the simulation using the mesh size h = 4 μm. It can be seen from these figures that the twins propagate along the same paths, regardless of the mesh size used in each simulation. During the twin propagation process, the thickness of the deformation twin evolves according to Eq. (73), which accounts for the temperature and the resolved shear stress on the twin system under investigation. Fig. 20 shows contour plots of the twin thickness at various points in time during the deformation process, where the field plotted in color is the thickness b of the embedded weak discontinuity contained within each element. It is noted that the patch of elements containing the embedded twin has serrated (zig-zag) boundaries due to the finite element spatial discretization. Also note that this patch of elements has a different thickness depending on the mesh size used. However, it can be seen from Fig. 20 that the average twin thickness is approximately 0.3 μm, irrespective of the (significantly larger) mesh sizes used in the two cases considered (2 μm and 4 μm). In other words, by comparing the results obtained from simulations using two different mesh sizes, we conclude that the evolution of the embedded twin thickness does not depend on the underlying mesh size. It is also emphasized that, consistent with experimental observations, the embedded twin follows a straight propagation path which intersects (i.e. passes through) the set of finite elements forming the aforementioned serrated patch. Fig. 21 compares the distributions of the σ 22 stress in the twinned and untwinned crystal regions after the deformation twin is fully developed. Numerical results obtained from simulations using different mesh sizes are compared, as done previously. It is clear from Fig. 21, that the middle grain has a relatively smooth stress distribution inside and outside the embedded twin, and that there are no visible stress concentrations around the patch of elements containing the embedded twin. 4. Summary In this paper, we present a new computational approach to incorporate deformation twinning as a plastic deformation mechanism, in addition to crystallographic slip, into an explicit finite element formulation. Different deformation gradients are adopted for the untwinned and twinned crystal regions based on the underlying kinematics, which is a departure from the conventional isodeformation gradient assumption that is widely used in the literature. More specifically, in the
24
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Fig. 20. Evolution of the twin thickness b within the elements that contain the embedded twin. Numerical results obtained using the mesh size h = 2 μm (a-c) are compared to those obtained using the mesh size h = 4 μm (d–f). Note that it is the value of b within each element, not the size h of the element itself, that determines the thickness of the embedded twin.
untwinned crystal region, plastic deformation is dominated by plastic slip and the deformation gradient is multiplicatively decomposed into elastic and plastic parts. This is consistent with large deformation crystal plasticity formulations. In the twinned crystal region, deformation twinning is considered as an additional plastic deformation mechanism, and the associated deformation gradient is multiplicatively decomposed into elastic, plastic and twinning parts. The two deformation gradients are related through the compatibility condition and the traction continuity condition at the interface between the untwinned and twinned regions. In the finite element formulation, the deformation twin is modeled as a heterogeneous deformation mode, and treated as weak discontinuities embedded within individual elements. Multiple twins are allowed to nucleate at the grain boundaries within a poly-crystalline structure by using a probability-based twin nucleation model. Once initiated at the grain boundaries, twin propagation and growth are simulated according to the corresponding evolution equations under subsequent deformation. Several numerical examples are provided to verify the implementation, and to demonstrate the capabilities of the proposed crystal plasticity finite element framework with regard to the accurate, mesh-independent representation of the deformation twinning process. The computational framework presented in this paper is only the first step toward incorporating deformation twinning in general three-dimensional (3D) finite element simulations of crystal plasticity. Future improvements can be achieved at least through the following avenues. The single-crystal elasto-viscoplasticity model adopted in this paper is based on a purely mechanical theory and assumes isothermal conditions, i.e. although the temperature is treated as a state variable, it is not allowed to evolve during the deformation process. In order to depict the underlying physics related with deformation twinning more accurately, a constitutive model that adequately represents thermo-mechanical effects will be adopted in future work. The current constitutive model can easily be replaced by such a coupled thermo-mechanical model without significant changes to the structure of the proposed computational framework itself. For all the simulations presented in this work, deformation twinning is simplified as a two-dimensional deformation process via spatial rotations of the crystal lattices under investigation. However, this plastic deformation mechanism, like crystallographic slip, is intrinsically threedimensional. Therefore, more sophisticated numerical methods, especially with regard to the evolution of the twin topology, are required to ensure the propagation continuity of the deformation twin in a general three-dimensional setup. Addition-
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
25
Fig. 21. Distributions of the stress component σ 22 in (a,b) the untwinned crystal regions and (c,d) the twinned crystal regions after the twin is fully developed in the middle grain (t = 60.0 × 10−9 s) using different mesh sizes.
ally, the propagation and growth of deformation twins are highly complex multi-physical processes that involve different length-scales. Combined theoretical and experimental research efforts are required to arrive at an accurate description of this physical phenomenon. Acknowledgments Funding for this work was provided by the Laboratory Directed Research and Development (LDRD) program of Los Alamos National Laboratory under project number 20150431ER, by the DoD/DOE Joint Munitions Program (JMP), and by the Advanced Simulation and Computing–Physics and Engineering Models (ASC–PEM) program. The authors wish to acknowledge fruitful discussions about the present work with Hansohl Cho, Biao Feng, Marko Knezevic, Veronica Livescu, and Jason R. Mayeur. Appendix A. Newton scheme for the traction continuity equation The traction continuity condition (37) can be recast as:
n · [[σ ]] = 0,
(A.1)
where [[σ [[ represents the jump of the stress field at the interface between the embedded twin and the surrounding material. Recall that during the time interval [tn , tn+1 ], t = tn+1 − tn represents the time increment, (•)n represents the variable
26
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
at time tn that is already known, and (• )n+1 represents the counterpart at time tn+1 that needs to be updated. The above nonlinear equation is to be solved using a Newton scheme. Let (•)(k) and (• )(k+1 ) represent the quantities associated with the k-th and (k + 1 )-th iterations. Based on Eq. (43), we have
δ FUn+1 = δ LU FUn ,
δ FZn+1 = δ LZ FZn .
(A.2)
By using the mid-point rule and the relationship between α Z and α U shown in Eq. (35), we recast the rate form between L and L, as shown in Eq. (30), into the incremental form:
LZi j = δik δ jl + (h/b − 1 )α U Ti j Tkl Lkl in the embedded twin, L i j = LUi j = δik δ jl − α U Ti j Tkl Lkl in the untwinned region.
(A.3)
ˆ , defined as follows, Introducing the variables αˆ and m
αˆ = α U Tkl Lkl , mˆ = αˆ m,
(A.4)
we rewrite Eq. (A.3) in component form:
LUi j = Li j − mˆ i n j ,
LZi j = Li j +
or in tensorial form:
LU = L − mˆ ⊗ n,
L Z = L +
h ˆ in j , −1 m b
(A.5)
h ˆ ⊗ n. −1 m b
Substituting into Eq. (A.2), we obtain
(A.6)
δ FUn+1 = δ LU FUn = −δ mˆ ⊗ n FUn , δ FZn+1 = δ LZ FZn =
h −1 b
δ mˆ ⊗ n FZn .
(A.7)
The Cauchy stress is linearized, yielding the following expression for the stress at the end of the (k + 1 )-th iteration at t = tn+1 : +1 ) ) ) ) σ n(k+1 ≈ σ n(k+1 + δσ n(k+1 = σ n(k+1 +
) ∂ σ n(k+1 ) ) : δ Fn+1 = σ n(k+1 + W n(k+1 : δ Fn+1 , ∂ Fn+1
(A.8)
where W represents the material Jacobian associated with the single-crystal elasto-viscoplasticity material model of Section 2.3. The derivation of the Jacobian W is presented in Appendix B. Using the linearization of the stress shown in Eq. (A.8), the traction continuity condition (A.1) is linearized as follows: +1 ) n · [[σ ]]n(k+1 =n ·
σZ − σU
(k+1) n+1
≈n ·
σ Z + δσ Z
(k )
n+1
−
σ U + δσ U
(k ) n+1
= 0.
(A.9)
For notational convenience, we drop the subscript indicating the time step and rearrange the above equation into the following form,
n·
δσ Z − δσ U
(k )
= −n ·
σZ − σU
(k )
(A.10)
According to Eqs. (A.7) and (A.8), we have
δσ Z = W Z : δ FZ = W Z : and
h −1 b
δ mˆ ⊗ n FZn ,
(A.11)
δσ U = W U : δ FU = −W U : δ mˆ ⊗ n FUn .
(A.12)
Substituting Eqs. (A.11) and (A.12) into (A.10), we obtain (in component form):
ni
h Z U − 1 WiaZ jk Fn,pk + WiaU jk Fn,pk b
( k )
ˆ j = −ni n pδm
σiaZ − σiaU
(k )
.
(A.13)
This can be recast into the following compact form:
ˆ j = g(ak ) . za(kj ) δ m
(A.14) (k )
where the (k + 1 )-th iteration residual ga
g(ak ) = −ni za(kj ) = ni
U (k )
σiaZ − σia
(k )
and tangent matrix za j are given by
,
h Z U − 1 WiaZ jk Fn,pk + WiaU jk Fn,pk b
(A.15)
( k ) n p.
(A.16)
ˆ . Finally, we recover the values of α U and The linearized system of Eq. (A.14) is used to solve iteratively for the vector m ˆ with the aid of Eq. (A.4). m from the converged solution for m
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
27
Appendix B. Derivation of the material Jacobian Our goal here is to derive the Jacobian of the single-crystal elasto-viscoplasticity model, defined as follows:
W (τ ) =
∂ σ (τ ) . ∂ F (τ )
This material Jacobian is used to solve the nonlinear traction continuity condition, using the Newton scheme described in Appendix A. Quantities at the beginning of the time step (•)t are known. Dropping the time indicator τ associated with quantities at the end of the time step for convenience, and expressing all tensors in indicial notation, we have:
∂σi j ∂ 1 e e e = F S F ∂ Fkl ∂ Fkl detFe im mn jn e 1 e ∂ Fjne ∂ Fim e e 1 e ∂ Smn e e e e e e ∂ = S F + F F + F S + F S F , im im mn jn detFe ∂ Fkl mn jn ∂ Fkl jn im mn ∂ Fkl ∂ Fkl detFe
Wi jkl =
and
1 ∂ 1 ∂ Fe e −T = − F : . ∂ Fkl detFe detFe ∂ Fkl
(B.1)
(B.2)
Letting
T =
∂ Fe ( τ ) , ∂ F (τ )
Q=
∂ Se ( τ ) , ∂ F (τ )
(B.3)
we have
Wi jkl =
1 e e e e e e e e e e −1 Timkl Smn Fjn + Fim Qmnkl Fjn + Fim Smn T jnkl − Fim Smn Fjn Fqp T pqkl . e detF
(B.4)
• Calculation of T :
Ti jkl =
∂ Fiej ∂ p = Fim Fmn (t )−1 δn j − γ α Sα0n j ∂ Fkl ∂ Fkl α
∂ Fim p −1 p α α −1 ∂ α α = F (t ) γ S0n j + Fim Fmn (t ) γ S0n j δn j − δ − ∂ Fkl mn ∂ Fkl n j α α
∂ γ α p p −1 α α −1 α δn j − = δik δml Fmn (t ) γ S0n j + Fim Fmn (t ) − S ∂ Fkl 0n j α α
p = δik Flnp (t )−1 δn j − γ α Sα0n j − Fim Fmn (t )−1 Rαkl Sα0n j , α
(B.5)
α
where
Rαkl =
∂ γ α . ∂ Fkl
• Calculation of Q:
Se (τ ) = Se trial (τ ) −
(B.6) α
γ α C α ( τ ) ,
Qi jkl =
∂Cα ∂ Siej ∂ Siejtrial α ∂ γ α = − Ci j − γ α i j = Di jkl − Ciαj Rαkl − γ α Jiαjkl , ∂ Fkl ∂ Fkl ∂ F ∂ F kl kl α α α α
Di jkl =
∂ Siejtrial , ∂ Fkl
(B.7)
(B.8)
where
Jiαjkl =
∂ Ciαj . ∂ Fkl
(B.9)
Let A = Ce trial (τ ) = Fe trial (τ )T Fe trial (τ ), thus,
Se trial = C : Ee trial = C :
Di jkl =
1 e trial 1 C − 1 = C : [A − 1], 2 2
1 ∂ Siejtrial ∂ 1 ∂ Amn 1 = Ci jmn (Amn − δmn ) = Ci jmn = Ci jmn Lmnkl , ∂ Fkl ∂ Fkl 2 2 ∂ Fkl 2
(B.10)
(B.11)
28
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
where
Lmnkl =
∂ Amn . ∂ Fkl
Recall
(B.12)
Cα = C : sym Ce trial Sα0 = C : sym[ASα0 ] =
1 T C : ASα0 + Sα0 A , 2
(B.13)
therefore,
∂ Ciαj 1 ∂ 1 = C A Sα + Sα0 pm A pn = Ci jmn Lmpkl Sα0 pn + Sα0 pm L pnkl . ∂ Fkl 2 ∂ Fkl i jmn mp 0 pn 2
Jiαjkl =
(B.14)
• Calculation of L: T
Fe trial (τ ) = F(τ )Fp (t )−1 ,
Ce trial (τ ) = Fe trial (τ ) Fe trial (τ ) = Fp (t )−T F(τ )T F(τ )Fp (t )−1 ,
(B.15)
∂ Ai j ∂ p −1 ∂ p p = Fmi (t ) Fnm Fnp Fppj (t )−1 = Fmi (t )−1 [F F ]F (t )−1 ∂ Fkl ∂ Fkl ∂ Fkl nm np p j p = Fmi (t )−1 δnk δml Fnp + Fnm δnk δ pl Fppj (t )−1
Li jkl =
p = Flip (t )−1 Fkm Fmp j (t )−1 + Fmi (t )−1 Fkm Flpj (t )−1
(B.16)
• Calculation of Rα :
γ α = γ˙ α (τ α , sαρ ) t = γ α (τ α , sαρ ), (B.17) where the slip rate γ˙ α is a function of the resolved stress τ α and the slip resistance sα . For the sake of simplicity, following ρ the assumption that the change in the slip resistance sα ρ during the time step is neglected (Balasubramanian, 1998; Dai, 1997), we have γ α = γ α (τ α ). For metallic materials, the elastic stretch is usually infinitesimal, therefore,
τ α = (Ce Se ) : Sα0 ≈ Se : Sα0 = τ α (Se ) ⇒ γ α = γ α (Se ),
(B.18)
and
Rαi j =
∂ γ α ∂ γ α ∂ Skle = = Bαkl Qkli j , ∂ Fi j ∂ Skle ∂ Fi j
Bαkl =
∂ γ α . ∂ Skle
(B.19)
Substituting the above equations into the calculation of Q, expressed in Eq. (B.8), we get:
Qi jkl = Di jkl − = Di jkl −
α
α
Ciαj Rαkl −
α
γ α Jiαjkl
Ciαj Bαmn Qmnkl −
or in tensorial form:
Q=D−
( Cα ⊗ Bα ) : Q −
α
α
α
γ α Jiαjkl ,
(B.20)
γ α J α .
(B.21)
Rearranging the terms in Eq. (B.21), we can solve for Q as follows:
Q=K
−1
:
D−
• Calculation of Bα :
Bα =
α
γ α J α
,
Ki jkl = δik δ jl +
α
Ciαj Bαkl .
(B.22)
∂ γ α ∂ γ α ∂τ α ∂ γ α ∂ e α ∂ γ α = = S : S = sym[Sα0 ]. ( ) 0 ∂ Se ∂τ α ∂ Se ∂τ α ∂ Se ∂τ α
(B.23)
Based on the derivations from Eq. (B.1) to Eq. (B.23), and recalling that and F = F(t + t ) are known, the computational procedure for calculating W (τ ), defined in Eq. (65), is summarized as follows: Fe (t),
p
p
p
p
1. Li jkl = Fli (t )−1 Fkm Fm j (t )−1 + Fmi (t )−1 Fkm Fl j (t )−1 , 2. Di jkl = 12 Ci jmn Lmnkl , 3. Jiαjkl = 12 Ci jmn [Lmpkl Sα + Sα L ], 0 pn 0 pm pnkl ∂ γ α
4. Bα = ∂τ α sym[Sα ], ij 0i j 5. Ki jkl = δik δ jl + α Ciαj Bα , kl α α 6. Qi jkl = Ki−1 ( D − mnkl α γ Jmnkl ), jmn αQ , 7. Rα = B ij kl kli j p p 8. Ti jkl = δik Fln (t )−1 (δn j − α γ α Sα ) − Fim Fmn (t )−1 α Rαkl Sα0n j , 0n j 9. Wi jkl =
1 [T Se F e detFe imkl mn jn
e Q e e e e e e e −1 + Fim mnkl Fjn + Fim Smn T jnkl − Fim Smn Fjn Fqp T pqkl ].
Fp (t),
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
29
References Acharya, A., Beaudoin, A., 20 0 0. Grain-size effect in viscoplastic polycrystals at moderate strains. J. Mech. Phys. Solids 48 (10), 2213–2230. Addessio, F.L., Luscher, D.J., Cawkwell, M.J., Ramos, K.J., 2017. A single-crystal model for the high-strain rate deformation of cyclotrimethylene trinitramine including phase transformations and plastic slip. J. Appl. Phys. 121 (18), 185902. Agrawal, V., Dayal, K., 2015. A dynamic phase-field model for structural transformations and twinning: regularized interfaces with transparent prescription of complex kinetics and nucleation. part i: formulation and one-dimensional characterization. J. Mech. Phys. Solids 85, 270–290. Agrawal, V., Dayal, K., 2015. A dynamic phase-field model for structural transformations and twinning: regularized interfaces with transparent prescription of complex kinetics and nucleation. part II: two-dimensional characterization and boundary kinetics. J. Mech. Phys. Solids 85, 291–307. Anand, L., 2004. Single-crystal elasto-viscoplasticity: application to texture evolution in polycrystalline metals at large strains. Comput. Methods Appl. Mech. Eng. 193 (48), 5359–5383. Ardeljan, M., Beyerlein, I., McWilliams, B., Knezevic, M., 2016. Strain rate and temperature sensitive multi-level crystal plasticity model for large plastic deformation behavior: application to AZ31 magnesium alloy. Int. J. Plast. 83, 90–109. Ardeljan, M., McCabe, R., Beyerlein, I., Knezevic, M., 2015. Explicit incorporation of deformation twins into crystal plasticity finite element models. Comput. Methods Appl. Mech. Eng. 295, 396–413. Asaro, R., Needleman, A., 1985. Texture development and strain hardening in rate dependent polycrystals. Acta Metall. 33, 923–953. Asaro, R., Rice, J., 1977. Strain localization in ductile single crystals. J. Mech. Phys. Solids 25 (5), 309–338. Asaro, R.J., 1983. Micromechanics of crystals and polycrystals. In: Advances in Applied Mechanics, 23. Elsevier, pp. 1–115. Balasubramanian, S., 1998. Polycrystalline Plasticity: Application to Deformation Processing of Lightweight Metals. Massachusetts Institute of Technology. Belytschko, T., Fish, J., Engelmann, B.E., 1988. A finite element with embedded localization zones. Comput. Methods Appl. Mech. Eng. 70, 59–89. Belytschko, T., Liu, W.K., Moran, B., Elkhodary, K., 2013. Nonlinear finite elements for continua and structures.. John wiley & sons. Beyerlein, I., Capolungo, L., Marshall, P., McCabe, R., Tomé, C., 2010. Statistical analysis of deformation twinning in magnesium. Phil. Mag. 90, 2161–2190. Beyerlein, I., Tomé, C., 2010. A probabilistic twin nucleation model for HCP polycrystalline metals. Proc. Royal Soc. A 466, 2517–2544. Beyerlein, I.J., Tomé, C.N., 2008. A dislocation-based constitutive law for pure Zr including temperature effects. Int. J. Plast. 24, 867–895. Bronkhorst, C., Hansen, B., Cerreta, E., Bingert, J., 2007. Modeling the microstructural evolution of metallic polycrystalline materials under localization conditions. J. Mech. Phys. Solids 55 (11), 2351–2383. Busso, E.P., 1990. Cyclic deformation of monocrystalline nickel aluminide and high temperature coatings. Massachusetts Institute of Technology. Capolungo, L., Marshall, P., McCabe, R., Beyerlein, I., Tomé, C., 2009. Nucleation and growth of twins in Zr: a statistical study. Acta Mat. 57, 6047–6056. Cheng, J., Ghosh, S., 2017. Crystal plasticity finite element modeling of discrete twin evolution in polycrystalline magnesium. J. Mech. Phys. Solids 99, 512–538. Cheng, J., Shen, J., Mishra, R.K., 2018. Discrete twin evolution in Mg alloys using a novel crystal plasticity finite element model. Acta Mat. 149, 142–153. Chester, S.A., Bernier, J.V., Barton, N.R., Balogh, L., Clausen, B., Edmiston, J.K., 2016. Direct numerical simulation of deformation twinning in polycrystals. Acta Mat. 120, 348–363. Clayton, J.D., 2011. Nonlinear Mechanics of Crystals. Springer. Clayton, J.D., Knap, J., 2011. A phase field model of deformation twinning: nonlinear theory and numerical simulations. Physica D 240, 841–858. Clayton, J.D., Knap, J., 2016. Phase field modeling and simulation of coupled fracture and twinning in single crystals and polycrystals. Comput. Methods Appl. Mech. Eng. 312, 447–467. Dai, H., 1997. Geometrically-necessary dislocation density in continuum plasticity theory, FEM implementation and applications. Massachusetts Institute of Technology. Dequiedt, J.L., Denoual, C., Madec, R., 2015. Heterogeneous deformation in ductile FCC single crystals in biaxial stretching: the influence of slip system interactions. J. Mech. Phys. Solids 83, 301–318. Feng, B., Bronkhorst, C., Addessio, F., Morrow, B., Cerreta, E., Lookman, T., Lebensohn, R., Low, T., 2018. Coupled elasticity, plastic slip, and twinning in single crystal titanium loaded by split-hopkinson pressure bar. J. Mech. Phys. Solids 119, 274–297. Fish, J., Belytschko, T., 1988. Elements with embedded localization zones for large deformation problems. Comp. Struct. 30, 247–256. Ghosh, S., Cheng, J., 2018. Adaptive multi-time-domain subcycling for crystal plasticity FE modeling of discrete twin evolution. Comput. Mech. 61, 33–54. Gray III, G.T., 1997. Influence of strain rate and temperature on the structure-property behavior of high-purity titanium. J. Phys. IV France 7. C3–423–C3–428 Gurao, N.P., Kapoor, R., Suwas, S., 2011. Deformation behavior of commercially pure titanium at extreme strain rates. Acta Mat. 59 (9), 3431–3446. Gurtin, M.E., Fried, E., Anand, L., 2010. The Mechanics and Thermodynamics of Continua. Cambridge University Press. Hansen, B.L., Bronkhorst, C.A., Ortiz, M., 2010. Dislocation subgrain structures and modeling the plastic hardening of metallic single crystals. Modelling Simul. Mater. Sci. Eng. 18, 055001. Hill, R., Rice, J., 1972. Constitutive analysis of elastic-plastic crystals at arbitrary strain. J. Mech. Phys. Solids 20 (6), 401–413. Hughes, T.J.R., 20 0 0. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover Publications. Jin, T., Mourad, H.M., Bronkhorst, C.A., 2019. A comparative study of shear band tracking strategies in three-dimensional finite elements with embedded weak discontinuities. Finite Elem. Anal. Des. 155, 11–31. Jin, T., Mourad, H.M., Bronkhorst, C.A., Livescu, V., 2018. Finite element formulation with embedded weak discontinuities for strain localization under dynamic conditions. Comput. Mech. 61, 3–18. Jin, T., Mourad, H.M., Bronkhorst, C.A., Livescu, V., Zhang, X., Linder, C., Regueiro, R.A., 2019. Three-dimensional explicit finite element formulation for shear localization with global tracking of embedded weak discontinuities. Comput. Methods Appl. Mech. Eng. 353, 416–447. Kalidindi, S., Bronkhorst, C., Anand, L., 1992. Crystallographic texture evolution in bulk deformation processing of FCC metals. J. Mech. Phys. Solids 40 (3), 537–569. Kalidindi, S.R., 1998. Incorporation of deformation twinning in crystal plasticity models. J. Mech. Phys. Solids 46 (2), 267–290. Knezevic, M., Lebensohn, R., Cazacu, O., Revil-Baudard, B., Proust, G., Vogel, S.C., Nixon, M., 2013. Modeling bending of α -titanium with embedded polycrystal plasticity in implicit finite elements. Mater. Sci. Eng. A 564, 116–126. Kocks, U., 1976. Laws for work-hardening and low-temperature creep. J. Eng. Mater. Technol. 98 (1), 76–85. Kothari, M., Anand, L., 1998. Elasto-viscoplastic constitutive equations for polycrystalline metals: application to tantalum. J. Mech. Phys. Solids 46, 51–83. Kröner, E., Teodosiu, C., 1974. Lattice defect approach to plasticity and viscoplasticity. In: Sawczuk, A. (Ed.), Problems of plasticity. Noordhoff International Publishing, Leiden, pp. 45–82. Kubin, L.P., 2013. Dislocations, mesoscale simulations and plastic flow. Oxford University Press. Kumar, M.A., Beyerlein, I.J., McCabe, R.J., Tomé, C.N., 2016. Grain neighbor effects on twin transmission in hexagonal close-packed materials. Nature Commun. 7, 13826. Kumar, M.A., Beyerlein, I.J., Tomé, C.N., 2016. Effect of local stress fields on twin characteristics in HCP metals. Acta Mat. 116, 143–154. Kumar, M.A., Capolungo, L., McCabe, R.J., Tomé, C.N., 2019. Characterizing the role of adjoining twins at grain boundaries in hexagonal close packed materials. Sci. Rep. 9, 3846. Kumar, M.A., Kanjarla, A.K., Niezgoda, S.R., Lebensohn, R.A., Tomé, C.N., 2015. Numerical study of the stress state of a deformation twin in magnesium. Acta Mat. 84, 349–358. Leclercq, S., Nguy, C., Bensussan, P., 1989. Microstructure transformations in pure polycrystalline α -titanium under static and dynamic loading. Mechanical Properties of Materials at High Rates of Strain. Institute of Physics, Bristol. Lee, E., 1969. Elastic-plastic deformation at finite strains. J. Appl. Mech. 36 (1), 1–6. Lee, E.H., Liu, D.T., 1967. Finite-strain elastic–plastic theory with application to plane-wave analysis. J. Appl. Phys. 38 (1), 19–27.
30
T. Jin, H.M. Mourad and C.A. Bronkhorst et al. / Journal of the Mechanics and Physics of Solids 133 (2019) 103723
Levitas, V.I., Preston, D.L., 2005. Thermomechanical lattice instability and phase field theory of martensitic phase transformations, twinning and dislocations at large strains. Phys. Lett. A 343, 32–39. Liu, C., Shanthraj, P., Diehl, M., Roters, F., Dong, S., Dong, J., Ding, W., Raabe, D., 2018. An integrated crystal plasticity-phase field model for spatially resolved twin nucleation, propagation, and growth in hexagonal materials. Int. J. Plast. 106, 203–227. Livescu, V., Beyerlein, I.J., Bronkhorst, C.A., Dippo, O.F., Ndefru, B.G., Capolungo, L., Mourad, H.M., 2019. Microstructure insensitive twinning: a statistical analysis of incipient twins in high-purity titanium. Materialia 6, 100303. Luan, Q.M., Britton, T.B., Jun, T.S., 2018. Strain rate sensitivity in commercial pure titanium: the competition between slip and deformation twinning. Mater. Sci. Eng. A 734, 385–397. Mourad, H.M., Bronkhorst, C.A., Livescu, V., Plohr, J.N., Cerreta, E.K., 2017. Modeling and simulation framework for dynamic strain localization in elasto-viscoplastic metallic materials subject to large deformations. Int. J. Plast. 88, 1–26. Niezgoda, S.R., Kanjarla, A.K., Beyerlein, I.J., Tomé, C.N., 2014. Stochastic modeling of twin nucleation in polycrystals: an application in hexagonal close– packed metals. Int. J. Plast. 56, 119–138. Peirce, D., Asaro, R., Needleman, A., 1982. An analysis of nonuniform and localized deformation in ductile single crystals. Acta Metall. 30 (6), 1087–1119. Rapperport, E., Hartley, C., 1960. Deformation modes of zirconium at 77 K, 300 K, 575 K, and 1075 K. Trans. Am. Inst. Min. Metall. Eng. 218, 869. Rice, J., 1971. Inelastic constitutive relations for solids: an internal-variable theory and its application to metal plasticity. J. Mech. Phys. Solids 19 (6), 433–455. Salem, A.A., Kalidindi, S.R., Doherty, R.D., 2003. Strain hardening of titanium: role of deformation twinning. Acta Mat. 51 (14), 4225–4237. Salem, A.A., Kalidindi, S.R., Doherty, R.D., Semiatin, S.L., 2006. Strain hardening due to deformation twinning in α -titanium: mechanisms. Metall. Mater. Trans. A 37 (1), 259–268. Simo, J.C., Hughes, T.J.R., 1986. On the variational foundations of assumed strain methods. J. Appl. Mech. 53, 51–54. Song, S., III, G.G., 1995. Influence of temperature and strain rate on slip and twinning behavior of Zr. Metall. Mater. Trans. A 26 (10), 2665–2675. Tjahjanto, D.D., Turteltaub, S., Suiker, A.S.J., 2008. Crystallographically based model for transformation-induced plasticity in multiphase carbon steels. Continuum Mech. Thermodyn. 19 (7), 399–422. Van Houtte, P., 1978. Simulation of the rolling and shear texture of brass by the Taylor theory adapted for mechanical twinning. Acta Metall. 26 (4), 591–604. Wang, J., Hirth, J., Tomé, C., 2009. (1¯ 012 ) twinning nucleation mechanisms in hexagonal-close-packed crystals. Acta Mat. 57, 5521–5530. Wang, J., Hoagland, R., Hirth, J., Capolungo, L., Beyerlein, I., Tomé, C., 2009. Nucleation of a (1¯ 012 ) twin in hexagonal close-packed crystals. Scripta Mater. 61, 903–906. Wang, J., Beyerlein, I.J, Tomé, C., 2010. An atomic and probabilistic perspective on twin nucleation in Mg. Scripta Mater. 63, 741–746. Wang, J., Yadav, S., Hirth, J., Tomé, C., 2013. Pure-shuffle nucleation of deformation twins in hexagonal-close-packed metals. Mater. Res. Lett. 1, 126–132.