A finite element methodology to incorporate kinematic activation of discrete deformation twins in a crystal plasticity framework

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 358 (2020) 112653 www.elsevier.com/locate/cma

A finite element methodology to incorporate kinematic activation of discrete deformation twins in a crystal plasticity framework Matthew Kasemera ,∗, Paul Dawsonb b

a Department of Mechanical Engineering, The University of Alabama, Tuscaloosa, AL 35487, United States of America Sibley School of Mechanical and Aerospace Engineering, Cornell University, Ithaca, NY 14850, United States of America

Received 31 March 2019; received in revised form 22 September 2019; accepted 23 September 2019 Available online xxxx

Abstract A finite element framework is presented in which grains within a polycrystalline microstructure are pre-discretized into lamellar regions that become candidates to deform by twinning. This enables the mapping of the motion due to twinning on geometrically proper regions. At some point in the loading of the polycrystal, a twin is triggered, at which point the nodes within a twin region are rapidly mapped to their twinned location, the region’s crystal lattice is reoriented, and the remainder of the body deforms by means of crystallographic slip to enforce mechanical equilibrium. The framework allows for both the onset of twinning and twin growth, and is intended to serve as a test bed for investigating models proposed for these phenomena. In this paper, the framework is described and simulations are performed to demonstrate its ability to handle the large, fast deformation of twinning within a targeted grain of a polycrystalline aggregate. Various simulations are performed, where the initial twin width and the point of twin activation in the loading history are parameterized. Deformation fields in and around the twinned region are inspected after a twinning event, and changes in local stress states are discussed in light of global and local energetic metrics. c 2019 Elsevier B.V. All rights reserved. ⃝ Keywords: Crystal plasticity; Finite element method; Deformation twinning; Discrete twinning

1. Introduction Metals and metallic alloys are chosen for many engineering applications with demanding structural requirements. Understanding their deformation responses is critical for the effective use of the alloys in their design applications and for successfully manufacturing components fabricated with them [1]. Metals, however, are often complex materials with complex behaviors, and understanding their deformation responses is non-trivial [2,3]. Primarily, heterogeneity in terms of a metal’s microstructure, as well as complex deformation response at the crystal level [4,5], contribute to limited understanding of both its elastic and plastic responses. In light of these characteristics, predicting deformation responses of metals quickly becomes unmanageable without the aid of mathematical models [6,7]. Simulating the deformation of materials by representing their microstructures provides insight into the development of plasticity at the crystal scale [8–10], and the influence that both the microstructure and the crystal ∗ Corresponding author.

E-mail address: [email protected] (M. Kasemer). https://doi.org/10.1016/j.cma.2019.112653 c 2019 Elsevier B.V. All rights reserved. 0045-7825/⃝

2

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

behavior have on the resulting macroscopic behavior [11]. Understanding the microstructure–property relationship of polycrystalline metals is best achieved through the use of sophisticated crystal scale finite element simulations performed on high-fidelity representations of microstructures. Metallic single crystals may appear in many different atomic configurations, chiefly crystals that exhibit either cubic or hexagonal symmetry. Hexagonal crystals are especially interesting due to their presence in high strength structural alloys (such as titanium alloys [1]), as well as their relatively high degree of anisotropy in both their elastic and plastic responses [12–19]. Over a broad range of temperature and strain rate, both cubic and hexagonal metals deform plastically primarily by means of crystallographic slip. Slip is a deformation mode in which shearing results from the relative displacements of atomic planes within a crystal, and is accomplished by means of the net motion of dislocations along crystallographic planes and directions. Slip deformation tends to be diffuse (spread over an entire grain), but exhibit heterogeneous patterns [10]. Slip propagates at a rate determined by the local stress, and can cause shear of an arbitrary magnitude. Lending to even further complexity, crystals may exhibit other modes of plastic deformation in addition to crystallographic slip [2,4,5]. Deformation twinning is a different mode of plastic deformation which presents as characteristically distinct from slip in nearly every regard [20,21]. The shear associated with twinning is defined and uniform over a discrete volume, rather than arbitrary and diffuse. Additionally, the lattice within a twin experiences a large, fixed reorientation, rather than evolving incrementally depending on the local state. Although much is known regarding the geometry of twins, the understanding of the kinetics of twinning lags behind other mechanisms like slip [2,11,22–25]. It has been observed that twins develop at very high shear rates compared to slip [2,22,24,25], but rate (kinetic) equations are uncertain. Further, the condition(s) needed to activate a twin are not well-established, and remain an active topic of research. Upon activation, the appearance of a twin has the potential to quickly change local stress states in the vicinity of the twin, as large deformation occurs to accommodate the twin, coupled with the introduction of new grain boundaries (twin boundaries) [11,26]. The introduction of new boundaries and reorientation of the crystal lattice could have effects upon the strength and ductility of the aggregate [21]. While various crystal types may exhibit twinning under certain conditions, the discussion herein will focus entirely on hexagonal crystals due to their propensity to twin. The dissimilarities between slip and twinning have been an impediment to including twinning within simulations of the deformations of polycrystalline metals, as noted in a review of modeling efforts devoted to polycrystalline solids [11]. The discrete motion of twins is difficult to represent, as formulations that are designed for the continuous development of strain under slip often are not well-suited for a large fixed-size jump in strain from twinning. Additionally, the merging of two deformation modes, twinning and crystallographic slip, with potentially very disparate time scales due to their shear rate disparity [2,22,24,25], renders a system analogous to a stiff differential equation [27], lending to further computational complexity. This resulted in the development of models which trade away the prediction of local responses to instead focus on accurate prediction of bulk responses [28–30]. Various models have been proposed to predict the deformation response due to twinning through judicious relaxation of one or more of these physical attributes of twinning. Chief among them, models which consider twinning as a “pseudo-slip” system have emerged as the primary method to represent twinning as a deformation mode [30–32]. Twin systems are considered at each material point as additional slip systems. The model homogenizes the deformation responses due to both slip and twinning at a material point by considering the “volume fractions” of the material assumed to be deforming by either deformation mode (that is, a homogenization variable determining the contribution of each deformation mode toward the total deformation response). This scheme does not fully consider the physical differences between the two deformation modes, as the discrete nature of twinning is not incorporated, and the disparity between the shear rates associated with slip and twinning is obscured. Thus, the prediction of local stress fields due to twinning is hampered. The ability to discretize a grain into multiple elements – and thus the abandonment of homogenization techniques in these types of simulations – has led to this class of model being employed beyond its intended scope. Computational capabilities have increased such that the resolution of polycrystal deformation simulations now include fine microstructural features [33,34]. Grains are now commonly discretized into many elements, which allows for complex local heterogeneous deformation to be studied [10,35]. Instantiation methods have similarly been developed to rapidly and reliably construct realistic polycrystal geometries and accompanying meshes [33,36]. Further, crystal scale finite element frameworks [37,38] have advanced to be both fast and reliable — as parallelization of code has allowed for much larger simulations to be performed, and confidence in

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

3

existing models has grown with respect to their ability to model plastic deformation by means of crystallographic slip [39]. Indeed, the majority of contemporary simulations of the deformation of crystalline solids are performed using crystal plasticity models on three-dimensional representations of microstructures [11] — a departure from simulations performed even as late as the mid 2000s, where models were often simplified [40], and microstructural representations were either low fidelity [8] or reduced to relatively small material volumes [41]. With these advances, current frameworks are poised to approach the problem of deformation twinning in ways previously hindered. Owing to these recent developments, other twinning modeling efforts have been implemented in a wide variety of frameworks, including visco-plastic self consistent (VPSC) frameworks [42–44], which accounts for some neighborhood effects, though does not explicitly model the full-field microstructure. More recently, attempts have been made at including discrete deformation twinning in full-field frameworks. Namely, twinning models have been incorporated into fast Fourier transform (FFT) [45–47] and finite element modeling (FEM) [48–52] frameworks. Among finite element frameworks considering discrete deformation twinning, the work of Cheng et al. [51,52] and Ardeljan et al. [48–50] are especially notable. Cheng et al. focus on the proposal of physics-based models for twin nucleation and growth — perhaps the most complex twinning model to date. The work of Ardeljan et al. concerns attempts to consider discrete twins via lamellar geometries, though at the cost of adaptive re-meshing and mapping at the point of twin activation. Their results indicate interesting trends concerning stress distributions — specifically stress increases in neighboring grains in the vicinity of the twin tip. It is apparent from the continuing struggle to successfully integrate twinning into the computational approaches now widely used for slip-dominated plasticity that progress has been hindered by the fact that twinning presents several inter-related challenges stemming from its unique attributes. Separating effects of the geometry of twinning from the nucleation of a twin, and that from the kinetics of the twinning process, is difficult. However, without being able to separate these factors in model development, proposed models for the separate portions of the twinning process (i.e., nucleation and growth) cannot be readily assessed. The present paper offers a modeling framework that can be used as a research platform to study the consequences of assumptions on the predicted responses during twinning events. This is accomplished by taking advantage of a principal attribute of the finite element method — namely the interpolation functions associated with the mapping of the motion that lie at the core of its implementation. Having local control over the motion facilitates enforcing the precise geometry of deformation twinning within a grain. With this accomplished, the consequences of the twinning deformation on the stress fields within and surrounding the twin can be computed and critiqued. In this way, the framework temporarily decouples critical aspects of twinning, allowing one to be fixed while the other varied for the purpose of model assessment. We expect that the framework will serve as a flexible test bed for studying twin nucleation and growth models (the framework serves as an open platform for model implementation), and potentially for studying displacive transformation induced plasticity (e.g., martensitic transformation) where similar modeling challenges prevail [53]. To demonstrate the framework capabilities, a single grain at the center of a three-dimensional polycrystalline aggregate is pre-discretized via a multilevel tessellation method [35] into intragrain lamellar regions which may be considered for twinning. Lamellar regions are oriented to represent a single twin system. For the sake of demonstration, only a single twin system per grain is considered — though the framework is poised to handle multiple variants per grain. Multiple instantiations are generated with various initial twin thicknesses. A finite element mesh comprised of well-formed elements is generated for the resulting geometry, allowing for the consideration of elements and nodes associated with an appropriately shaped twin lamella. Twins are activated at pre-selected points in the deformation path, one in the elastic domain and one beyond the initial yield point. No specific twin nucleation models are considered in this demonstrative example. Motion due to twinning is achieved by mapping the nodal points contained within a discrete twin region to their expected twinned location at a relatively rapid rate. Lattice reorientation is updated by rotating lattice orientations within all elements of the twin volume to their expected twinned orientations. The finite element framework determines the response of the surrounding volume, allowing it to deform as needed via crystallographic slip to satisfy the field equations. Crucial aspects of the mechanical response of the aggregate are computed, including stress distributions and energy metrics such as the total and plastic work. Results are viewed as a function of initial twin width and the point of twin activation in the loading history. Specifically, the changes to stress fields and the energy associated with twinning events are closely inspected. It is envisioned that the new framework will be used in conjunction with closely coordinated experiments. In particular, recent developments in experimental capabilities present researchers with possible paths to study

4

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

Fig. 1. (a) Illustration of a single crystal, where the shaded area indicates some volume deformed by twinning, and the circle highlights the large lattice reorientation due to twinning (i.e., lattice mirrors across the twin plane, dashed line). (b) Illustration of a single crystal, where the shaded area indicates a region undergoing a similar shear via slip, and the circle highlights a slip boundary. Reorientation is relatively minor compared to that associated with twinning, and is consequently not illustrated. Table 1 Primary slip systems of hexagonal close packed crystals. Slip system name

Number of systems

m

s

Basal Prismatic Pyramidal

3 3 12

{0 0 0 1} {0 0 1 0} {0 0 1 1}

⟨1 1 2 0⟩ ⟨1 1 2 0⟩ ⟨1 1 2 3⟩

deformation twinning more deeply [54,55]. High energy X-ray diffraction (HEXD) has emerged as a potentially viable technique to observe the formation of twins during in situ loading. The development of fast digital detectors has allowed for the precise measurement of diffraction peaks associated with different grains within a microstructure [9], and the ability to track crystals which may be twins associated with other crystals [54,56]. The recent advent of high energy diffraction microscopy, or “near field” diffraction, allows for a non-destructive mapping of a material’s microstructure [55,57,58]. In conjunction with one another, and expected advancement in both techniques, experimental observation and understanding of deformation twinning are likely to accelerate [55,56,59]. 2. Kinematics 2.1. Plastic deformation of crystals Plastic deformation is an isochoric process. In general, six strains are necessary to describe a strain state, though the isochoric constraint reduces the number of independent strains to five. Plastic deformation can be accomplished by a number of physical mechanisms. At low and moderate temperatures, slip and twinning often are the most active mechanisms [22]. For plastic deformation to be capable of achieving any shape change by means of slip, a total of five linearly independent slip systems are necessary [2]. In crystals exhibiting hexagonal symmetry, three families of slip systems are often considered to meet this requirement. The most common families considered are the basal, prismatic, and pyramidal slip systems. The basal and prismatic slip systems are referred to as ⟨a⟩-type slip systems — as they allow the crystals to slip only in the direction of one of the ⟨a⟩ axes. Similarly, pyramidal slip systems are referred to as ⟨c + a⟩-type slip systems, as they allow for slip in a direction with a component aligned with the ⟨c⟩-axis. These slip systems are summarized in Table 1. Each slip family has an associated slip strength — that is, the shear stress sufficient to sustain plastic shearing associated with a slip system (the crystallographic combination of a slip plane and slip direction) by means of the net movement of dislocations [60]. Slip is commonly observed on ⟨a⟩-type slip systems due to their ease of activation (low strength). In some hexagonal close packed crystals, the ⟨c + a⟩-type systems may be difficult to activate due to a relatively high strength [1]. This relegates slip to directions perpendicular to the ⟨c⟩-axis, thus inhibiting general plastic deformation if only slip is considered. In these metals, deformation twinning may be observed.

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

5

Table 2 Description of {1 0 1 2} twin system of hexagonal close packed crystals. K1

η1

K2

η2

S

{1 0 1 2}

⟨1 0 1 1⟩

{1 0 1 2}

⟨1 0 1 1⟩

0.167

Fig. 1 illustrates a single crystal specimen sheared due to twinning, juxtaposed next to a single crystal specimen similarly sheared but due to slip. Note the difference between both deformation modes. During twinning, a volume of material undergoes a uniform shear of fixed magnitude (the “twin”, shaded in gray), whereas the rest of the material remains undeformed (the “parent”). Additionally, the crystal lattice experiences a large reorientation in the twin such that it mirrors across the twin plane (depicted as a dashed line). In contrast, during slip a dislocation shifts one plane of atoms relative to another. The accumulated motion of many dislocations creates arbitrary amount of strain that may be heterogeneously spread over the deforming volume. The amount and direction of lattice reorientation depends on the particular slip systems activated. Additionally, the shear rate due to twinning is potentially much faster (several orders of magnitude) than that associated with slip [2,11,22,24,25]. Twin systems allow for the necessary deformation in the direction of the crystal ⟨c⟩-axis [1]. A primary example is the {1 0 1 2} twin system, though there are other twin systems which are witnessed to activate in hexagonal crystals [1]. Twin systems may be broadly described by the direction of twin shear, η1 (described by a direction ⟨h k i l⟩), and the plane on which shearing occurs, K 1 (described by a plane {h k i l}). In addition to the direction of twin shear and the twin plane, twins may also be described in terms of a second undeformed plane (K 2 in the configuration pre-twinning, K 2′ in the configuration post-twinning), a direction perpendicular to the first undeformed plane on the second undeformed plane (η2 and η2′ in similar configurations), their width, w, and the shear magnitude, S. The geometry of the {1 0 1 2} twin system is summarized in Table 2, and a schematic of the geometry may be viewed in [1,61]. The reorientation of the lattice is described as a 180◦ rotation of the crystal lattice about the twin plane normal — or, in other words, a reflection across the twin plane [21,28]. For hexagonal titanium crystals, twinning on the {1 0 1 2} system produces a reorientation such that the ⟨c⟩-axes of the parent crystal and twinned crystal deviate by ∼85◦ [1]. For example, a parent crystal with its ⟨c⟩-axis aligned with the loading direction, upon twinning, reorients such that its ⟨c⟩-axis is nearly perpendicular to the loading direction. This substantial reorientation of the crystal lattice has the potential to markedly affect the local response of the material. The directional elastic response of hexagonal crystals is strongly dependent on the orientation of the crystal ⟨c⟩-axis with respect to the loading axis, and can exhibit a relatively large degree of anisotropy [2,12,39]. The large degree of plastic anisotropy with respect to crystallographic slip will also be significantly affected by such a large reorientation. 2.2. Kinematics of motion In general, the deformation of a material point is described via a deformation gradient composed of an elastic stretch, a rotation, and plastic deformation. The total deformation gradient for a material point is multiplicatively decomposed into these three components (elastic, Fe , rotation, F R , and plastic, F P ): F = Fe F R F P .

(1) e

The elastic portion of deformation is assumed to be composed of small elastic strains, E , and is approximated by: Fe = I + Ee where ∥Ee ∥ ≪ 1.

(2)

Classically, the plastic portion of deformation is considered to have a contribution only from the deformation associated with crystallographic slip. However, for materials exhibiting both slip and twinning, the plastic portion of deformation must be further decomposed into a portion considering crystallographic slip, Fsli p , and one considering deformation twinning, Ftwin , such that: F P = Fsli p Ftwin .

(3)

6

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

The deformation gradient for crystallographic slip develops from the shearing associated with active slip systems as discussed in the preceding section (and further described in Section 3.1). The deformation gradient associated with a deformation twin is defined as [61]: ( ) Ftwin = I + S dtwin ⊗ ntwin (4) where S is the shear magnitude, dtwin is the shear direction, and ntwin is the twin plane normal. Note that this describes a simple shear mode. The shear direction and the twin plane normal in Eq. (4) are Cartesian vectors. The shear direction is associated with the crystallographic direction, η1 , while the twin plane normal is associated with the crystallographic plane K 1 (as described in the preceding section and summarized in Table 2). The shear magnitude is constant for a given twin system. 2.3. Motion of a discrete region due to twinning The aim of this framework is to construct the ability to activate discrete deformation twins in a polycrystalline aggregate at relatively high shear rates. This is achieved by pre-discretizing a polycrystalline aggregate (illustrated in two dimensions in Fig. 2(a)) via a multilevel tessellation approach (described in detail in [35]) to include lamellar regions at the intragrain scale (Fig. 2(b)). Within a single grain, lamellar regions which may be considered for twinning are otherwise indistinguishable from other regions within the grain. A finite element mesh is constructed for the geometry including lamellar features (Fig. 2(c)). The identification and selection of elements and nodes within a lamellar region is thus facilitated via this meshing (Fig. 2(d)). At the onset of twinning, velocities are imposed on the selected region’s nodes (Fig. 2(e)), with velocities and time steps chosen to provide shear rates several orders of magnitude faster than those associated with crystallographic slip. The motion due to twinning is imposed on the twin region (Fig. 2(f)), and slip is allowed to occur throughout the aggregate, such that mechanical equilibrium is enforced. During the imposition of motion, the crystal lattice is reoriented within the region. Note that this approach is in contrast to the work of Ardeljan [48], in that the twin geometry is considered prior to simulation — allowing this framework to forego mid-simulation re-meshing at the onset of twinning. The specific sequence of steps to impose the motion due to twinning on a region is summarized in the following sections. Within a discrete region experiencing twinning, dtwin and ntwin are calculated for every element. As some misorientation will likely develop within a region prior to twinning, elements within a region will have different values for dtwin and ntwin owing to the spread in orientations, and thus different values for the calculated Ftwin . To ensure compatibility of elements within a discrete twin region, the same deformation gradient, ⟨Ftwin ⟩, is applied to all elements within a region: ( ) ⟨Ftwin ⟩ = I + S ⟨dtwin ⟩ ⊗ ⟨ntwin ⟩ (5) where ⟨dtwin ⟩ is the weighted average of the elemental twinning directions, and ⟨ntwin ⟩ is the weighted average of the elemental twin plane normals. This assumption is inconsequential if the development of misorientation within a region is small — valid at lower strains, certainly at strains below macroscopic yield. For all nodes belonging to elements within a twin region, the nodal displacement, utwin , is calculated as the difference of the twinned nodal locations as calculated via the average deformation gradient due to twinning and the pre-twinned nodal locations: ( ) utwin = ⟨Ftwin ⟩x − x (6) where x are the current nodal coordinates prior to the onset of twinning. To be clear, this ensures that the deformation prior to a twinning event is considered. The total displacement is achieved incrementally over a number of steps, n steps : uincr =

utwin n steps

(7)

where each step has an equal time step, ∆t. The displacement is applied incrementally to both aid in convergence of the solver, and avoid very large stress jumps as the motion is imposed. Finally, nodal velocities, vtwin , are calculated: vtwin =

uincr . ∆t

(8)

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

7

Fig. 2. Illustrations depicting a two-dimensional sample (a) discretized into grain structures, with (b) further discretization of a single grain into intragrain lamellar structures, and (c) an attendant finite element mesh of the geometry including lamellar structures. (d) Identification of mesh features (elements and nodes) within a twin region, (e) velocities imposed on nodes, and (f) nodes of twin region mapped to expected twinned locations. Lamellar regions may not be separated by strictly parallel lines due to tessellation regularization [33].

( ) Time steps are chosen to be relatively small, such that the shear(rate associated with the twinning event, O 10 s−1 , ) is relatively rapid in comparison to the nominal strain rate, O 0.01 s−1 — though this is adaptable for instances where a twinning event occurs much slower. 3. Constitutive equations and state evolution It would now be customary to introduce the constitutive equations relating stress and deformation for each deformation mechanism — elasticity, plasticity by slip, and plasticity by twinning. However, as the framework is designed to kinematically impose the motion due to twinning, we negate the need for a stress–twinning relationship to facilitate solving the boundary value problem. Instead, we impose the twinning motion on the twin’s volume, and since the motion is completely imposed on the twinning volume to precisely map a twinning event, it is possible to prescribe the stress within the twinning volume rather than having to solve for it through the complete finite element methodology. This is the same for any problem for which the motion is fully prescribed: stresses are evaluated from the constitutive equations for a known motion. The material surrounding the twinned region is forced to conform at the twin interface. Stress in the surrounding material will be subject to the field equations, including constitutive equations for elasticity and slip, described in the following section. Thus, for the implementation of the framework, all of the body uses the same formulation except in the twins while they are twinning. In those volumes the stress is prescribed from the known motion. The base framework herein employs a rate dependent, restricted slip model [40] and Voce type hardening law [3] to predict the deformation of single crystals due to slip. A parallelized finite element framework is utilized as the solution method. A description of the finite element implementation is summarized in [62]. The model described in this section contains no inherent physical length scale, though it is able to capture relative grain size effects [35] (note: the propensity for a grain to twin may be correlated to the size of the grain [63], though these effects are not explicitly considered in this study). Various studies have shown the ability of the framework to accurately predict both bulk response such as texture evolution [3], and local response such as complex intragrain heterogeneous

8

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

deformation [10,11]. Importantly, results suggest that the spatial arrangement of grains and phases have large effects on the development of the aggregate deformation response both locally and globally, and must be accurately represented to understand the plastic response at the crystal scale [8,10]. 3.1. Constitutive equations Plasticity by restricted slip is governed by a rate-dependent model. To better accommodate this rate dependence in the numerical implementation, kinematics may also be considered in a rate form. To this end, the velocity gradient, L, is deduced from the deformation gradient: ˙ −1 L = FF

(9)

where F˙ denotes the time rate of change of the deformation gradient. The velocity gradient is decomposed into the deformation rate, D, and spin, W: L = D + W.

(10)

Substituting Eq. (1) into Eq. (9), and separating the deformation rate into its hydrostatic and deviatoric components yields expanded forms of the deviatoric deformation rate and the spin (note that the small elastic strain assumption in Eq. (2) has been employed): ′ ˆ P ′ + Ee′ W ˆ P −W ˆ P Ee′ D′ = E˙ e + D ˆ P + Ee′ D ˆ P′ − D ˆ P ′ Ee′ and W = W

(11)

where “′ ” denotes the deviatoric component of a given tensor. The plastic portions of the deviatoric deformation rate and the spin are transformed to the plastic intermediate configuration by means of the rotation portion of the deformation gradient: ˆ P′ = FR DP′ FR T D ˆ P = F R W P F R T + F˙ R F R T . and W

(12)

Constitutive equations for the elastic and plastic responses are now substituted into Eq. (11) to introduce the stress. The elastic portion of deformation is governed by Hooke’s law, τ = C(r)Ee

(13)

where the Kirchhoff stress, τ , is related to the elastic strain through the use of a fourth order anisotropic stiffness tensor, C, which is reduced to reflect the symmetry of the crystal [5]. Additionally, orientation dependence of the stiffness tensor is included by means of the Rodrigues parameterization, r, of a crystal orientation [64]. The plastic portion of the velocity gradient is decomposed into the plastic deformation rate and the plastic spin ˆ P ′ , is a linear in a similar fashion to Eq. (10). For plasticity by crystallographic slip, the plastic deformation rate, D combination of simple shearing modes: ∑ ˆ P′ = D γ˙ k Pˆ k (14) k ( ) ˆk . where Pˆ k = sym sˆk ⊗ m Here, γ˙ k denotes the shearing rate on a slip system, k, and Pˆ k is the symmetric portion of a slip system’s Schmid ˆ k . Slip systems commonly tensor — constructed using a slip system’s slip direction, sˆk , and plane normal, m considered for hexagonal close packed crystals are detailed in Table 1. The plastic spin is defined as: ∑ ˆk ˆ P = F˙ R F R T + W γ˙ k Q (15) k ( ) ˆ k = skw sˆk ⊗ m ˆk where Q ˆ k is the skew portion of the Schmid tensor. where Q

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

9

The shearing rate of a given slip system is defined using a power law expression to introduce rate dependence (controlled by m), and relates the shearing rate on a slip system to the resolved shear stress of a system, τ k : ( k ) m1 |τ | k γ˙ = γ˙0 sgn(τ k ) (16) gk k k ′ where τ = tr(Pˆ τ ). Here, the shearing rate on a slip system is scaled by the fixed-state strain rate scaling coefficient, γ˙0 , and is dependent on a slip system’s current strength, g k . In this study, we assume m to be constant and equal across slip families, though recent studies suggest more complex behavior in hexagonal crystals [65,66]. 3.2. State evolution The state of a given point in a material is described using the current crystallographic orientation and the current slip system strengths. The evolution equations of these variables are cast in a rate form. First, the evolution of the orientation of a crystal is again written in terms of its Rodrigues parameterization: 1 (ω + (ω · r) r + ω × r) 2 where ω is the lattice spin, defined as: ( ) ∑ P k k ˆ . ˆ − ω = vect W γ˙ Q r˙ =

(17)

(18)

k

Second, evolution of each slip system strength is defined using a Voce hardening assumption: ( k )n gs (γ˙ ) − g k ∑ k g˙ = h 0 h kl |γ˙l | gsk (γ˙ ) − g0k k ( ∑ ⏐ ⏐ )m ′ ⏐ k⏐ k γ˙ k . where gsk (γ˙ ) = gs0 γ˙s0

(19)

Here, the matrix, h kl , controls the hardening interaction for slip system “k” with other slip systems “l”. The instantaneous rate of evolution is controlled by the strength hardening rate coefficient, h 0 , the current slip system strength, g k , the initial strength, g0k , the saturation strength, gsk , and the non-linear Voce exponent, n. Furthermore, the saturation strength depends on deformation rate according to a power-law relation with the net rate of slip, γ˙ . Some studies indicate that the onset of twinning modifies the hardening behavior inside the twin region [11,67,68], though these effects are ignored in the present study. 4. Sample generation and twin activation 4.1. Discretization of twin regions The cornerstone of the extended framework presented in this study is the ability to discretize single grains into regions which may be considered as twin regions at some point in deformation, as summarized at the beginning of Section 2.3 (see: Fig. 2). Upon discretization of a domain into grains and possible twin regions, an attendant finite element mesh is constructed. The deformation due to twinning may be applied, thus, to discrete lamellar regions, maintaining a key difference between crystallographic slip and twinning. In short, a multilevel tessellation method is employed to discretize a domain into a number of grains and to further discretize grains into regions which may be considered for deformation twinning. A finite element mesh is constructed for the resulting geometry [33]. Various experimental studies of deformation twins [4,55,69] have witnessed twins exhibiting a lamellar geometric morphology. While twins of different shapes are known to form (e.g., lenticular twin regions), this present study focuses solely on deformation twins of the lamellar type. Note that lenticular twins could potentially be considered via the geometric discretization described in this section, by mapping motion not to an entire lamellar region, but by masking the geometry of a lamella to consider a lenticular geometry — to be clear, utilizing the same geometric

10

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

Fig. 3. A domain discretized into (a) grains, and (b) further discretization into candidate twin regions. Note that each grain is only discretized into twin regions representing only a single twin system, though each grain can be further discretized to represent other systems.

discretization. Focusing on lamellar twins allows for the use of existing frameworks to instantiate samples with twin regions. The multilevel tessellation method described in [35] established the ability to discretize domains into lamellae by means of the insertion of successive parallel planes into arbitrary domains. Finite element meshes are generated for these complex regions. Grains are discretized prior to simulation to include lamella which represent regions which may be considered for deformation twinning. Fig. 3 details an illustrative three-dimensional polycrystal with its grains discretized into lamellar twin regions, with each grain discretized considering only a single twin system (that is, only a single set of lamellae per grain). A first level of tessellation is generated to represent grains in a polycrystal (Fig. 3(a)). For this study, grain size and shape distributions are chosen for a Laguerre tessellation to generate equiaxed grains [34]. The second tessellation level is to generate lamellar regions within the grains of the first tessellation level to represent twin regions (Fig. 3(b)). Twin regions are generated using planar slices, which are defined by the twin width and the plane normals. In this illustrative example, each grain is discretized using only a single set of lamellar planes — representing only a single twin system. In general, however, to represent multiple twin systems, each grain would be discretized using multiple sets of lamellar planes. Within a grain, planar slices are roughly parallel, though may not be strictly so due to tessellation regularization [33] — necessary to ensure the generation of a robust mesh. Each grain’s crystallographic orientation must be known – as well as the twin system under consideration – to calculate the twin plane normal in the sample reference frame. Here, knowledge of the reciprocal lattice vector is employed to calculate the twin plane normal given the Miller–Bravais values for a given twin plane. The reciprocal basis, b, is first constructed from the crystal basis. With this basis, the reciprocal lattice vector is calculated with the Miller–Bravais indices for the crystallographic plane of interest, {h k i l}: ghkl = hb1 + kb2 + lb3 .

(20)

A plane normal in the original rectangular Cartesian basis is defined in the crystal frame as the normalized reciprocal lattice vector for that crystallographic plane: ghkl . (21) nChkl = ∥ghkl ∥ S This vector may be transformed to the sample coordinate system to arrive at nhkl . Procedurally, to instantiate lamellar twin regions within a polycrystalline aggregate, each grain is assigned a crystallographic orientation. The twin plane normals are then calculated for each discrete orientation. A tessellation is created, in which each grain, having been assigned a discrete orientation from the orientation set, is discretized

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

11

into twin regions using the calculated twin plane normals. Each lamellar region within a grain is assigned the same orientation. Thus, the lamellar features are not considered pre-existing microstructural features, as together they represent a single grain. A domain can be divided using any number of sets of planes described by their normals (that is, any number of twin variants). The increased subdivision of a grain as multiple twin systems are represented leads to the generation of successively finer geometric features. For a mesh of good quality (i.e., well formed elements), the inclusion of finer geometric features necessitates a higher element count. The increase in the number of calculations due to the higher element count is offset due to modern computational capabilities and expected further developments thereof. For the demonstrative simulations described later in this study, only a single grain will be discretized into a single twin system, though this does not represent the general case or suggest a limitation in the framework to only one twin system per grain. The twin system chosen for this present study is the (1 0 1 2)[1 0 1 1] twin system, and twin plane normals are calculated with this in mind. As a note, refining the twin lamellar thickness is essential in determining the response of the aggregate, in the same way that mesh refinement is essential in determining a converged solution. For instance, if lamellar regions are too thick, the point in deformation at which twinning is activated may be over-predicted (depending on the model used for activation), or the stress field over-predicted if a twin of unreasonably large thickness is activated too early. Decreasing the twin thickness, however, necessitates the ability to “grow”, or thicken the twin upon increased loading, as the stress field may be under-predicted upon continued loading if the conditions are proper for twin growth. The framework presented herein has the ability to consider twin growth via the activation of twinning in lamella neighboring a region which has already twinned. Activation of neighboring lamella would depend on the activation mechanism considered. Polycrystal generation, meshing, and visualization is achieved via the Neper software package [33,36]. 4.2. Imposition of twin motion and lattice reorientation To apply the motion due to twinning to a specific region, the nodal points within a region are considered. The nodal points of a discrete twin region are moved to their expected twin locations by temporarily imposing velocities on the nodes within a twin region such that the nodal points are properly mapped to their expected locations. The calculation for these velocities is described earlier in Section 2.3. For an ordinary simulation – such as the uniaxial tension of a domain in the shape of a rectangular prism discretized into a finite element mesh (Fig. 4(a)) – velocity boundary conditions are applied to two opposing surfaces. So called “grip” boundary conditions, for example, are intended to hold the area of the two opposing surfaces fixed, while one of these surfaces is held motionless and the other is the extension surface. Fig. 4(b) shows the applied grip boundary conditions on the surface nodes. The nodes on the bottom surface are fixed in all directions, while the nodes on the top surface have velocities applied such that the top surface extends away from the bottom surface yet is fixed in the other directions. For simulations without twinning, nodes which are given applied boundary conditions retain them throughout the course of the simulation. To apply the motion of a region due to deformation twinning, however, this restriction is temporarily relaxed. Velocities are imposed on all nodes belonging to a selected region — those that reside both within the volume, and on the boundary of a region to be twinned. Velocities are held constant on these nodes for a set period of time until the nodes within the region reach their expected twinned locations. During this time, all other nodes with applied boundary conditions are temporarily held in a fixed position (that is, all velocities at these nodes are temporarily set to zero). Note that the sample does not unload during this period, although there may be a slight relaxation of stress owing to the rate dependence of plastic flow if the stress is sufficiently high. This is, however, assumed to be small due to the short time period in which the boundary conditions are temporarily halted. Fig. 5(a) shows a single crystal discretized into lamellar regions, and its corresponding finite element mesh is shown in Fig. 5(b). A twin region of interest is highlighted in gray. Grip boundary conditions are initially applied to the mesh (Fig. 5(c)). At some point in deformation when a twin is activated, the grip conditions are temporarily held fixed, while velocities are imposed over any twin region (Fig. 5(d)). The pre-discretization of the domain into twin regions allows for the selection of nodes that belong to a lamellar region — or a region having the proper geometric morphology to represent a deformation twin. The mesh shown in Fig. 4(a) would not allow this, as it is formed with no knowledge of the orientation of the crystal or the twin system in question. When the twin motion is

12

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

Fig. 4. (a) Mesh of a single crystal instantiation, and (b) grip boundary conditions. Red dots represent fixed nodes and blue arrows represent applied velocities (arbitrary magnitude). For sake of clarity, boundary conditions are shown only at the element vertices.

complete, boundary conditions are returned to their state prior to the twinning event (Fig. 5(c)), and the simulation continues. The crystal lattice within a discrete twin region experiences a large reorientation. During a twinning event, the orientation of the twin is calculated based on the current orientation and the twin system in question. The reorientation of the lattice is described as a reflection of the crystal lattice about the twin plane — or, mathematically, a 180◦ rotation about the twin plane normal. As the Rodrigues parameterization for a rotation of 180◦ is indeterminate, it is necessary to consider the orientation of the crystals via a separate parameterization. Here, the reorientation of the crystal due to twinning is described via a rotation matrix [70]. This rotation is defined [28] via the direction cosines of the twin plane normal with the sample basis, ν, as: Rtwin = 2 (ν ⊗ ν) − I.

(22)

Orientations are considered for each element in the mesh, and thus each element within a twin region experiences reorientation due to twinning. The current orientation is updated to the twinned orientation, which includes both the rotation due to twinning and any lattice reorientation accrued prior to a twinning event. This reorientation is imposed on the first incremental twin step Eq. (7) and held fixed over the ensuing twin steps. This is chosen as other reorientation schemes – namely reorienting the lattice on the final twin step or incrementally reorienting the lattice – produce qualitatively and quantitatively comparable results. 4.3. Calculation of work The framework can inform the user of the work implications of the chosen model assumptions. If the consequences of certain model assumptions, say in an extreme case, violate the second law of thermodynamics, it is possible to reveal this via calculation and investigation of work. In general, the framework offers the determination of work from all of the mechanisms for evaluation by the user. Work is calculated both locally (on an element by element basis), or globally (over the entire domain of the polycrystal). The total work and the plastic work are calculated, and thus insight into the work due to elastic loading is inferred. Specifically, the work over the time steps during twin activation is investigated in an effort to gain an understanding of the amount of work necessary to activate a twin, and the effects felt locally in the vicinity of the twin. The work rate for a given deformation step of a simulation is defined using work conjugate pair [7] of the Cauchy stress and the deformation rate. Eq. (23) details the calculation of the (total) work rate for each element following

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

13

Fig. 5. (a) Single crystal discretized into lamella, with example twin region highlighted, (b) its mesh, (c) grip boundary conditions, and (d) temporary velocities imposed at the nodes during twin activation. Note how the velocities at the center-line of the twin lamella are near zero, where deformation would be minimal. Red dots represent fixed nodes and blue arrows represent applied velocities. For sake of clarity, velocities are shown only at the element vertices, and only on surface nodes (twin velocities are applied throughout the twin volume).

standard integration over an element by means of Gauss quadrature. W˙ =

n qp ∑

wi Ji (σ i : Di )

(23)

i

Here, n qp is the number of Gauss quadrature points in the element. Additionally wi is the quadrature weight, Ji is the determinant of the Jacobian, σ i is the Cauchy stress, and Di the deformation rate tensor, all evaluated at the quadrature points of a given element, i. The expression σ i : Di is the inner product between the Cauchy stress and the deformation rate tensor. The total work rate for the entire polycrystal is calculated as: W˙ tot =

n el ∑

W˙ j

(24)

j

where n el is the number of elements. The work rate is calculated using the stress and deformation rate output at the end of a deformation step. The work for a given step is approximated given the time step, ∆t, using the trapezoidal rule: ) ∆t ( ˙ (25) Wtot + W˙ totn Wstep = 2 where the quantity W˙ totn is the work rate at the end of the previous deformation step. To calculate the plastic work over a given step, Eq. (23) is modified to use the plastic deformation rate tensor, DiP . The same summation and integration described in Eqs. (24) and (25) apply. With the work rate and the work known for every time step and for every element, the global and local changes in temperature may be approximated. Assuming that the work is entirely converted to heat, the work is related to a change in temperature as: W ∆T = (26) ρC p V where W is the work for a region, ρ is the density of the material, C p is the specific heat, and V is the volume of a region. 5. Simulations The simulations performed for this study are designed to shed light on the mechanical response of a polycrystalline aggregate to the activation of a discrete deformation twin. That is, the polycrystalline aggregate used

14

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

Fig. 6. (a) A domain divided into grains, (b) a single grain at the center of the aggregate (shown at its location in the domain with area of interest highlighted), which is further discretized into twin regions.

in these simulations is designed specifically to view the development of the stress fields in and around the twin both prior to and after activation of a twin. These demonstrative simulations are not intended to mimic a specific observed twinning event, but rather to demonstrate the type and extent of data that may be extracted from the framework. As we show, the severity of the perturbation to the stress field depends on the timing of the twinning event relative to the onset of slip. The results show sensitivity to twin thickness and the point of twin activation in the loading history, which have implications for development and verification of twinning models. In addition, the example shows that a comprehensive picture emerges containing information of the local stress concentrations, the distribution of deformation and work around the twin, and the macroscopic stress–strain response — all important for critiquing a proposed model’s performance or defining expectations for a new model. The polycrystalline aggregate is generated by employing a Laguerre tessellation within the framework described in [34]. The tessellation is constructed in a cubic domain of initial edge length l = 1 mm (note: as the crystal plasticity finite element framework employed in this study contains no inherent length scale, all units herein are not absolute, but are provided to give the reader a relative sense of dimension). The polycrystal distributions for grain size (normalized equivalent diameter) and shape (sphericity) are chosen as delta functions of values 1 (by definition) and 0.85, respectively, in an effort to create an equiaxed microstructure (Fig. 6(a)). The sample is constructed with 100 grains, and a single grain near the center of the aggregate is chosen for discretization into lamellae to serve as possible twin regions (Fig. 6(b)). The number of grains is chosen such that the center grain is sufficiently far from the boundaries so as to reduce any significant effects due to boundary conditions. A sample with 100 grains and a delta distribution for the normalized equivalent diameter has an average equivalent grain diameter of approximately 0.27 mm. Multiple instantiations are produced in which the lamellar width for the center grain is variable. Fig. 7 shows the center grain of the polycrystal (Fig. 6(a)) discretized into lamellar regions of widths w = 0.0025 mm, w = 0.005 mm, and w = 0.01 mm. This allows for the study of the aggregate reaction to twin size. Twin sizes are chosen as fractions of the average grain diameter to qualitatively mimic twins witnessed in various studies (e.g., [71]). The center grain orientation is set to: ⎛ ⎞ 0.0 r = ⎝ 0.0 ⎠ (27) −0.2679 which corresponds to a rotation of 30◦ about the crystal ⟨c⟩-axis. The crystallographic orientation is chosen such that the ⟨c⟩-axis is aligned with the tensile direction, and the twin direction is in the sample x–z plane for the {1 0 1 2} twin system. The orientations of the other grains in the aggregate are chosen at random — that is, the aggregate is otherwise assumed to have a uniform orientation distribution. While twinning is largely suppressed in Ti–6Al–4V due to its high aluminum content [1], its model parameters are chosen for this study. The elastic [12] and plastic [13] parameters for Ti–6Al–4V are known with a high degree

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

15

Fig. 7. Region of interest of grain shown in Fig. 6(b) discretized into lamellar regions of width (a) w = 0.0025 mm, (b) w = 0.005 mm, and (c) w = 0.01 mm. Table 3 Single crystal elastic constants for the α phase of Ti–6Al–4V. C11 (GPa)

C12 (GPa)

C13 (GPa)

C44 (GPa)

169.7

88.7

61.7

42.5

Table 4 Initial slip system strengths for the α phase of Ti–6Al–4V. g0,b (MPa)

g0, p (MPa)

g0,π (MPa)

g0, p g0,b

390

468

663

1.2 : 1.7

:

g0,π g0,b

Table 5 Plasticity parameters used for the α phase of Ti–6Al–4V. h 0 (MPa)

gs0 (MPa)

m

m′

γ˙0 (s−1 )

γ˙s0 (s−1 )

190

530

0.01

0.01

1.0

5 × 1010

of confidence [39], and thus their use facilitate the inspection of stress fields due to twinning. Additionally, the modeling parameters for Ti–6Al–4V are representative of a large class of hexagonal crystals. Specifically, the character of the directional elastic modulus, and relative initial slip system strengths are sufficiently generic for metals exhibiting hexagonal symmetry. The high relative strength of the pyramidal system (nearly twice that of the basal system) would – in similar alloys – facilitate the development of deformation twinning. Additionally, Ti–6Al–4V exhibits little hardening, thus removing the complexities associated with the transient behavior of slip strengths. Indeed, previous work [39] has demonstrated that for Ti–6Al–4V the hardening assumption has little effect on slip initiation and the initial development of plasticity. Consequently, an isotropic hardening assumption is employed in this study. Simulation parameters are detailed in Tables 3, 4, and 5. With regard to alloys exhibiting hexagonal symmetry (specifically titanium alloys), the crystal plasticity finite element framework employed in this study (and variations thereof) has been a popular modeling choice. The framework has been successfully employed in a number of studies, and validated against a number of experimental techniques. Specifically, it has been shown to be adept at predicting the elastic response (as elucidated via HEXD) [12], the initiation of slip [39,72], and the onset and development of plasticity [13,73]. Additionally, this framework has been used in various computational studies on titanium alloys [35,74], and further studies into the behavior of macroscopic property convergence in titanium alloys [75]. Results of these studies both give confidence in the existing finite element framework’s ability to model crystallographic slip, and the employ of the constitutive parameters used in this study. Grip boundary conditions – as described in Section 4.2 – are employed to strain the aggregate. Samples are extended in their “z” directions with velocities of 0.01 mm s−1 . Twinning is activated at 0.5% macroscopic strain. The strain state at which twin is activated is chosen at an arbitrary point in the macroscopic elastic regime. Twin

16

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

Fig. 8. Representative stress–strain curve obtained from simulating the extension of the polycrystalline aggregate. The solid circle denotes the point at which twinning is activated, shortly before the onset of the elastic–plastic transition. Additional results for twin activation at 1% strain (after the initial yielding, denoted by the solid triangle) are presented in Section 6.2.

Fig. 9. (a) Aggregate with slice plane outlined, and (b) slice through the center of the aggregate with approximate region of interest outlined, both detailing the deviation of the crystal ⟨c⟩-axis from the loading direction, θ.

activation is applied over a number of steps, n steps = 10, such that the shear over each step is O (1%). Time steps of ∆t = 0.001 s are chosen, such that the twin forms over a total time of t = 0.01 s. During twinning, the strain path for many elements is altered due to rapid changes in the deformation field. Consequently, a number of relatively small time steps are taken directly after twin activation to aid in convergence. Table 6 outlines the time steps used in simulations with twin activation at 0.5% macroscopic strain. Fig. 8 shows a representative simulated stress–strain curve with twin activation at 0.5% macroscopic strain (denoted by a solid dot). 6. Effects of twinning on stress distributions Fig. 9 details the distribution of ⟨c⟩-axis orientations on the polycrystal (specifically the angle, θ , between each crystal ⟨c⟩-axis and the loading direction), and on a slice through the center of the aggregate, exposing the grain discretized for twinning and outlining a region of interest. Fig. 10 details θ within this region of interest prior to and after the activation of a twin for the simulations performed with various twin widths. Three-dimensional visualizations of the grain after twinning are shown in Fig. 11, where the deformation is exaggerated to accentuate

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

17

Table 6 Time steps used in simulations with twin activation at 0.5% macroscopic strain. Step numbers

∆t (s)

ϵmacr o (%)

1–5

0.100

0.500

6–10 11 12 13 14–27

0.001 0.005 0.040 0.050 0.100

0.505 0.510 0.550 0.600 2.000

Fig. 10. Deviation of the crystal ⟨c⟩-axis from the loading direction, θ , on the region of interest shown in Fig. 9(b) for the state (a) prior to twinning (with neighboring grains used in later analysis labeled “1” and “2”), and post twinning for simulations with twin width of (b) w = 0.0025 mm, (c) w = 0.005 mm, and (d) w = 0.01 mm.

the deformation due to twinning on a complex, three-dimensional domain. In all plots, the twin region presents itself as a band across the grain of interest, displaying the substantial reorientation of the crystal by highlighting the deviation of the twin’s ⟨c⟩-axis from the parent’s. Similarly, Fig. 12 details the von Mises stress throughout the polycrystal, and on a slice through the center of the aggregate. The stress field after the twinning event is shown on the interior slice of the polycrystal looking at the “x–z” plane in Fig. 13, and on the “y–z” plane in Fig. 14 for the simulations performed with various twin widths. This allows for the visualization of the complex stress states that form around the entire perimeter of the twin. Note how the effective stress experiences sharp increases around the perimeter of the twin region in neighboring grains, which dissipates as the distance from the twin region is increased.

18

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

Fig. 11. Three-dimensional visualizations of the grain of interest directly after twin activation plotted on the deformed mesh (exaggerated 5x), where dark regions represent the parent grain and light regions represent the twin (similar color scale to Fig. 10). Results from simulations with twin width (a) w = 0.0025 mm, (b) w = 0.005 mm, and (c) w = 0.01 mm.

Fig. 12. (a) Aggregate with slice plane outlined, and (b) slice through the center of the aggregate with approximate region of interest outlined, both detailing the effective (von Mises) stress.

Post-processing visualization of variables plotted on the finite element mesh is achieved via the Paraview software package [76]. 6.1. Influence of twin thickness The effective stress in the parent grain around the twin region decreases during the onset of twinning (compare Figs. 13(b), 13(c), and 13(d) to the pre-twin state in Fig. 13(a)). To quantify this stress drop, the average effective stress within the twin region and parent grain directly after twinning is calculated. Table 7 shows values for the average von Mises stress in various regions of the polycrystal for each twin width considered. The average stresses in neighboring grains show little change overall — indicating that, in sum, there is a near-equal amount of stress decrease in the polycrystal to accompany the large stress increases in the vicinity of the twin tips. The prediction of stress concentrations in neighboring grains supports some experimental observations [71], as well as results from similar computational studies [46,48–50]. Note how as the twin thickness increases, the stress field in the surrounding grains has larger and farther reaching changes (Fig. 13). As stated earlier, the effect of the twin thickness upon nucleation could have an effect on the over- or under-prediction of stress fields. Depending on the nucleation criteria or analysis, the thickest twin may be unrealistic to activate. As such, lamellar refinement studies – and the ability to grow twins – would be necessary. Twin growth would be considered through the

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

19

Fig. 13. Effective (von Mises) stress on the region of interest shown in Fig. 12(b) for the state (a) prior to twinning, and post twinning for simulations with twin width of (b) w = 0.0025 mm (with arrows approximating the direction of motion for the twin region), (c) w = 0.005 mm, and (d) w = 0.01 mm. Table 7 Average von Mises stress within the twin region, parent grain, total grain, and neighboring grains prior to and directly after twinning, calculated for simulations with twin activation at 0.5% macroscopic strain. Simulation

σv M Twin (MPa)

σv M Parent (MPa)

σv M Total grain (MPa)

σv M Neighboring grains (MPa)

Pre-twin

805

810

810

652

w = 0.0025 mm w = 0.005 mm w = 0.01 mm

518 457 467

720 647 581

717 642 576

652 651 652

successive activation of lamella neighboring the previously twinned region. In this way, a thicker twin may be eventually activated, but on a more realistic deformation path (depending on nucleation and growth mechanisms). Inspecting the stress components reveals more information concerning the stress field after the twinning event. This is best viewed by transforming the stress states to a coordinate system aligned with the characteristic twin directions. A coordinate frame, etwin is chosen, where e1twin = dtwin and e2twin = ntwin . Fig. 15 shows the normal component of stress in the direction of twinning (the “11” component in the transformed frame) after the twinning event for the simulations performed with various twin widths. Areas of tension and compression are evident at the twin tips extending into neighboring grains, due to the opposing motion on either side of the twin region. As twin width is increased, the volume effected by this motion is increased. Larger areas of neighboring grains are witnessed to have higher tensile or compressive stresses than the surrounding aggregate. Viewing the resolved shear stress on the twin system after the twin event allows for a comparison to similar two-dimensional studies. In a study utilizing a fast Fourier transform solution method [77], it is shown that the

20

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

Fig. 14. Effective (von Mises) stress on a slice perpendicular to that shown in Fig. 12 for the state (a) prior to twinning, and post twinning for simulations with twin width of (b) w = 0.0025 mm, (c) w = 0.005 mm, and (d) w = 0.01 mm.

resolved shear stress on the twin system experiences a drop within the twin region, and a sharp increase near the twin tips in surrounding grains immediately after the twinning event. This is reproduced here utilizing a full-field finite element solver. Fig. 16 shows the resolved shear stress along a trace through the center of the twin in the twin direction for simulations performed with various twin widths. Boundaries of the twin regions are roughly shown with dashed lines. In general — the resolved shear stress on the twin system decreases within the twin region after a twin event, but experiences a sharp increase in neighboring grains along this trace. The decrease in the resolved shear stress on the twin system is dependent on the twin width. In general, as twin width is increased, the resolved shear stress on the twin system experiences a decrease of greater magnitude (compare the post-twin states in Figs. 16(a) and 16(c)). Specifically, for results considering the thinnest twin in Fig. 16(a), note how the resolved shear stress on the twin system within the bounds of the twin region experiences a relatively low drop along this trace — in some places almost negligible. For the thickest twin in Fig. 16(c), the resolved shear stress does not approach the values from the pre-twin state, and thus experiences a greater decrease. In [46], it is assumed the driving mechanism for twin formation is the resolved shear stress on the twin plane (employing a pseudo-slip approach for twin nucleation). They postulate that a drop in the resolved shear stress within the twin region indicates an aversion to twin thickening. Under this hypothesis, results here indicate that the thinnest twin possesses the greatest ability to grow thicker, as there are regions in the twin with almost no drop in resolved shear stress. Additionally, the sharp increases in the maximum resolved shear stress on twin systems in neighboring grains further support results from similar studies [44,46]. This could lead to the activation – or transmission – of twinning to neighboring grains. Twin transmission will of course be dependent on the crystallographic orientations of neighboring grains, as predicted in a related study on modeling deformation twinning via a modified VPSC framework [78].

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

21

Fig. 15. Component of normal stress in the direction of twinning shear, for the state (a) prior to twinning, and post twinning for simulations with twin width of (b) w = 0.0025 mm (with arrows approximating the direction of motion for the twin region), (c) w = 0.005 mm, and (d) w = 0.01 mm.

The effects of mesh refinement in the parent grain and twin region is also considered. The base geometries for the various twin widths are re-meshed to increase the number of elements in the aggregate, and especially within the twin regions. Table 8 details the number of elements in the aggregate, parent grain, and twin region for both the base mesh and the refined mesh. The refinement level is chosen such that the refined mesh contains multiple elements across the width of the twin. In simulations performed with fine meshes, the effective stress decreases in the twin region and the parent grain at the onset of twinning by similar amounts to the simulations performed using the coarse mesh. Table 9 summarizes the average von Mises stress in various regions of the polycrystal both before and after the onset of twinning. The level of mesh refinement does not significantly alter the calculated average stresses in various regions throughout the polycrystal, nor did it cause the character of the stress fields to significantly deviate from the those plotted in Figs. 13 and 14, and for sake of brevity these redundant figures are omitted. 6.2. Twin activation in the elastic–plastic regime Additional simulations are performed in which twins are activated at 1.0% macroscopic strain. This state is chosen as it is within the elastic–plastic transition, to contrast against twinning in the elastic regime (refer to Fig. 8). Time steps are similar to those described in Table 6, with the primary difference being 10 initial steps prior to twin activation (rather than 5) — all other time steps are the same relative to the point at which the twin is activated. The same twin thicknesses, grain orientations, and microstructure morphology utilized above are again considered.

22

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

Fig. 16. Resolved shear stress for the twin system along a line through the center of the twin in the twin direction, for the state post twinning for simulations performed with twin width of (a) w = 0.0025 mm, (b) w = 0.005 mm, and (c) w = 0.01 mm.

Table 8 Number of elements contained in the total mesh, the parent grain, and the twin region for all simulations and two levels of mesh refinement. Domain subset

w (mm)

n elem (Coarse)

n elem (Fine)

Total

0.0025 0.005 0.01

107,453 52,916 38,219

917,931 647,608 350,418

Parent grain

0.0025 0.005 0.01

29,868 11,252 5,653

73,385 29,882 14,549

Twin region

0.0025 0.005 0.01

496 394 292

1,071 822 721

Fig. 17 details the von Mises stress prior to and post twinning for twin inserted at 1.0%. When twinning is activated at higher strain, the effect it has on stress fields in neighboring grains is not as appreciable as when twinning is activated at 0.5% macroscopic strain. Note how large stress concentrations near the twin tips apparent in Fig. 13 are not apparent in Fig. 17. This indicates that the activation of twinning at higher strains – and specifically here

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

23

Table 9 Average von Mises stress within the twin region, parent grain, total grain, and neighboring grains prior to and directly after twinning, calculated for simulations with twin activation at 0.5% macroscopic strain and performed with the refined mesh. Simulation

σv M Twin (MPa)

σv M Parent (MPa)

σv M Total grain (MPa)

σv M Neighboring grains (MPa)

Pre-twin

805

809

809

645

w = 0.0025 mm w = 0.005 mm w = 0.01 mm

517 452 464

720 648 584

718 642 579

651 642 641

Fig. 17. Effective (von Mises) stress on the region of interest shown in Fig. 12(b) for the state (a) prior to twinning, and post twinning for simulations with twin width of (b) w = 0.0025 mm (with arrows approximating the direction of motion for the twin region), (c) w = 0.005 mm, and (d) w = 0.01 mm. Results for simulations with twin activation at 1.0% macroscopic strain. Note the difference in the scale limits to those in Fig. 13, necessary to portray the features of the stress field at higher load.

above the elastic–plastic transition – leads to less pronounced effects in neighboring grains than when twinning is activated at lower (elastic) strain. This is likely because the single crystal yield surface constrains the stress magnitude. At 1.0% macroscopic strain, the stress state of each element is likely nearing or on the yield surface. As such, it is not possible to attain a large contrast such as shown in Fig. 13. The average stress in the parent grain and the twin again decrease with the activation of the twin. Average stresses in various subsets of the polycrystal are summarized in Table 10. This is seen visually in Fig. 17, where the activation of twinning has a pronounced effect on the local stress state in the parent grain. On average, the stresses in neighboring grains again see little change, though this is not due to an equal amount of stress increase and decrease in the polycrystal (as witnessed when the twin is activated in the elastic regime). Instead, this is

24

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

Table 10 Average von Mises stress within the twin region, parent grain, total grain, and neighboring grains prior to and directly after twinning, calculated for simulations with twin activation at 1.0% macroscopic strain. Simulation

σv M Twin (MPa)

σv M Parent (MPa)

σv M Total grain (MPa)

σv M Neighboring grains (MPa)

Pre-twin

1436

1445

1445

1088

w = 0.0025 mm w = 0.005 mm w = 0.01 mm

1020 875 761

1344 1262 1169

1340 1250 1150

1086 1076 1071

Table 11 Change in work and plastic work over the time steps of twin activation, calculated for simulations with no twin activation or twin activation at 0.5% macroscopic strain. Simulation

W (mJ)

W P (mJ)

No twin w = 0.0025 mm w = 0.005 mm w = 0.01 mm

0.0064 0.0178 0.0380 0.0689

∼0 0.0190 0.0395 0.0661

likely due to the fact that plasticity is already onset in neighboring grains as described above — leading to a less pronounced impact in the grains due to the onset of twinning. 7. Effects of twinning on global responses 7.1. Energetic observations Deformation twinning is inspected in light of the amount of work required to impose the deformation associated with twinning over a discrete region. The total and plastic work rate is calculated at the end of every deformation step, and the work is integrated through time by means of a trapezoidal integration scheme (see: Section 4.3). The entire body is considered, rather than local subsets, so as to understand the energetic effect twin activation has on the polycrystal as a whole. Of interest is the amount of work taken to activate a twin, and with this in mind, the change in work and plastic work is calculated between load steps 5 and 6 (for simulations with twin activation at 0.5% macroscopic strain). Table 11 summarizes results from simulations associated with twin activation at 0.5% macroscopic strain. Note that the plastic work is calculated to be negligible in the simulation without twinning (that is, negligible crystallographic slip), while the total work is non-zero. This indicates the body is deforming nearly entirely elastically. For simulations in which twins are activated at 0.5% macroscopic strain, the amount of work (both total and plastic) necessary to activate a twin increases as twin width increases. This is expected, as larger volumes of material are displaced, and larger portions of the polycrystal are affected by the onset of twinning. Work is dominated by the work due to plasticity — which considers the plastic deformation rate tensor, and thus the deformation due to crystallographic slip. For each twin width, the plastic work accounts for the majority of the work calculated over the polycrystal. This implies that during the twin event, the body is deforming largely by crystallographic slip so as to accommodate the deformation due to twinning. Note how in some instances (w = 0.0025 mm and w = 0.005 mm) the plastic work is greater than the total work, which indicates that the elastic work is negative over the time the twin is activated. In other words, more elastic unloading occurs during the onset of twinning than due to continued elastic loading over the entire polycrystal. For the thickest twin (w = 0.01 mm), the plastic work is less than the total work, indicating a net positive elastic work, or an increase in elastic loading over the volume of the polycrystal. Note how this trend correlates to the trend in the average stress drop in the twin region, summarized in Table 7. The thickest twin (w = 0.01 mm) does not experience as great a stress drop as the next thickest (w = 0.005 mm). The calculated total and plastic work over the onset of twinning is calculated as a function of mesh refinement. These values are summarized in Table 12. Again, the mesh refinement has little effect on the overall response of the twin region and parent grain during the onset of twinning. Compare these results to those from the simulations performed with coarse meshes, summarized in Table 11.

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

25

Table 12 Change in work and plastic work over the time steps of twin activation, calculated for simulations with no twin activation or twin activation at 0.5% macroscopic strain and performed with the refined mesh. Simulation

W (mJ)

W P (mJ)

No twin w = 0.0025 mm w = 0.005 mm w = 0.01 mm

0.0063 0.0177 0.0377 0.0683

∼0 0.0190 0.0394 0.0656

Table 13 Change in work and plastic work over the time steps of twin activation, calculated for simulations with no twin activation or twin activation at 1.0% macroscopic strain. Simulation

W (mJ)

W P (mJ)

No twin w = 0.0025 mm w = 0.005 mm w = 0.01 mm

0.0053 0.0114 0.0220 0.0407

0.0470 0.0678 0.0932 0.1285

The work over the course of the twinning event is similarly calculated for the simulations where twinning is activated at 1.0% macroscopic strain. These results are summarized for each simulation in Table 13. Note here that plasticity is dominant in every case, likely as the polycrystal is experiencing flow throughout most of the domain. It becomes less apparent at this stage whether the activation of the twin alone accounts for the large negative elastic work contribution. Compared against the values in Table 13, it is seen that – again, in this specific case – the activation of twinning requires more work than if the polycrystal were to deform entirely by slip, though this is not alarming for the above described reasons. Due to the nature of the integration technique, work calculations may be overestimated if abrupt elastic unloading occurs. A switch from a positive work rate to a negative work rate over a deformation step causes error in the integrated work. Numerical integration techniques will fail to capture the true negative work over a step, as the switch from a positive work rate to a negative work rate creates a saw-tooth curve, which is difficult to integrate using established techniques. This may be remedied by taking smaller time steps, thus finer resolving the work rate curve and allowing for more accurate integration. Consequently, the total work values presented in Tables 11 and 13 are likely overestimates. The calculation of work serves to provide a possible path to determine when a twin event may occur. The framework allows for the possibility of two simulations to be performed concurrently — one in which the body deforms at certain points by means of deformation twinning, and one in which the body deforms entirely by crystallographic slip. Integrated values for various metrics (e.g., work, elastic relaxation, or others) may be compared between the two simulations at each deformation step, and the mode of deformation that is energetically favorable may be chosen as the path which the simulation continues on. For the simulations performed, activation of a twin of any size requires more work than in simulations without a twinning event. Indeed, the total work is an order of magnitude higher over the twin step. This may be expected, as the twin is activated well before the onset of macroscopic plasticity. Localized plastic deformation has yet to become substantial at this point in loading, as the plastic work calculated for simulations without a twinning event is consistently negligible over that time step. A large deformation that relies on plasticity to accommodate the shape change may not compete in terms of an energetic measure to a body deforming nearly entirely elastically. Additionally, in these simulations, a twin is activated in an arbitrary crystal, embedded in an arbitrary aggregate, at an arbitrary point in deformation. No care is taken to simulate a “real” event. As such, it cannot necessarily be expected that the work calculations would indicate twinning as the favorable deformation mode (as opposed to crystallographic slip). This, however, provides a map for possible indication of twin nucleation. For instance, we see that the when a twin is activated in the elastic regime, the total work is equal to the plastic work, indicating little to no strain energy dissipation, compared to the activation of twinning in the elastic–plastic regime which shows a much larger strain energy dissipation, indicating that twinning is perhaps less favorable in the elastic regime.

26

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653 Table 14 Maximum local temperature change, volume of material experiencing maximum local temperature change, and average change in temperature for the entire polycrystal during twin activation, calculated for simulations with twin activation at 0.5% macroscopic strain. w (mm) 0.0025 0.005 0.01

Local ∆T (◦C)

V (mm3 )

Global ∆T (◦C)

43.7 50.5 36.4

9.7 × 10−7

0.0074 0.0159 0.0288

3.3 × 10−7 1.9 × 10−8

The change in temperature for a polycrystal at large is calculated per Eq. (26). Table 14 summarizes the average change in temperature over the entire polycrystal due to a twinning event, as well as the maximum local change and the volumes over which the maximum local changes occur. As twin width increases, the average temperature increase in the polycrystal increases. The maximum local increase in temperature is seen to be far greater than the average increase in temperature over the entire polycrystal. This is expected, as the twinning event – and the work associated with the event – is highly localized. Results indicate that the maximum temperature increase may be locally significant. This could necessitate the investigation of temperature dependent constitutive relations and evolution equations, though the magnitude of the maximum temperature increase is still relatively low, even in this worst-case scenario (assuming all work is translated to heat). 7.2. Slip activity in neighboring grains To study the effects that twin activation has on neighboring grains, the slip activity within the grains directly at the tip of the twin is investigated. The orientation of these two grains are such that their ⟨c⟩-axes deviate from the loading direction of the polycrystal by θ1 = 81.7◦ and θ2 = 26.25◦ , respectively — that is, one grain (grain 1) closely oriented to the twin region, and another (grain 2) more closely oriented to the parent grain. Slip activity is calculated as the percent of the grain volume containing active slip systems for a given family (basal, prismatic, or pyramidal). Fig. 18(a) details a summary of the slip activity in each grain for each simulation performed with twin activation at 0.5% macroscopic strain, and Fig. 18(b) details the same for each simulation performed with twin activation at 1.0% macroscopic strain. Slip activity is calculated at the load step directly after the twin event. Note that the sum of slip activity for a simulation may be greater than 100%, due to the possibility of polyslip in certain regions. When a twin is activated at 0.5% macroscopic strain, the volume in which slip activity is witnessed is an order of magnitude smaller than when a twin is activated at 1.0% macroscopic strain (observe the differences in the plot limits between Figs. 18(a) and 18(b)). This indicates that at 0.5% macroscopic strain – despite the large stress concentration in neighboring grains – plasticity is localized to only a small portion of these grains. Additionally, the amount of activity decreases as the twin width is decreased, which is consistent with the observations made concerning the stress fields. When a twin is activated at 1.0% macroscopic strain, the effects on slip activity in neighboring grains is less pronounced. In grain 1, for instance, the slip activity actually sees an overall decrease (perhaps due to some stress relaxation in portions of these grains) — though pyramidal slip systems are found to be active after twinning, while no pyramidal activity is observed when a twin is not activated. In grain 2, there again seems to be an increase in the amount of pyramidal activity as compared to the simulation in which no twin is activated. The slip families that show activity also differ between grains and activation point. Grain 1 exhibits nearly equal slip activity across all three slip families at 0.5% macroscopic strain, but almost entirely prismatic slip at 1.0% macroscopic strain. Conversely, grain 2 relies almost entirely on basal and pyramidal slip at both strains, but exhibits some (minor) prismatic activity at 0.5% macroscopic strain, while none is witnessed at 1.0% macroscopic strain. In other words, the slip activity in grains affected by twin activation is not only dependent on the orientation of the grain, but also the point at which twinning is activated.

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

27

Fig. 18. Slip activity directly after twin activation in the neighboring grains located at the twin tips (as labeled in Fig. 10(a)). Slip activity is presented as the percent volume of the grains deforming via each slip family, for simulations with twin activation at (a) 0.5% macroscopic strain, and (b) 1.0% macroscopic strain. Note the difference in scales between (a) and (b).

7.3. Macroscopic response The activation of a twin within the body has an appreciable effect on the macroscopic response of the aggregate. Fig. 19 details portions of the macroscopic stress–strain curves at the point at which twinning is activated for simulations performed with a twin width of w = 0.01 mm, and for twin activation at both 0.5% and 1.0% macroscopic strain. Note that in both cases, activation of twinning leads to a decrease in the tensile stress, and the stress–strain curve exhibits serrated behavior, expected from the formation of a twin [2,60]. The case of w = 0.01 mm is plotted, as it exhibits the most pronounced decrease in stress. Table 15 summarizes the stress decrease for all simulations, both in terms of magnitude, and as a percentage of the macroscopic tensile stress directly prior to twin activation. As the response is seen in the elastic regime (prior to flow in the aggregate at large) and beyond the plastic regime, the response is not due to relaxation during the momentary hold on strain during the twin event (as described in Section 4.2). Instead, the decrease in stress can be attributed to both the reorientation of the crystal lattice in the twin

28

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

Fig. 19. Macroscopic stress–strain curves for simulations performed with a twin width w = 0.01 mm, and twin activation at (a) 0.5% macroscopic strain, and (b) 1.0% macroscopic strain. In each plot, the inset details the stress–strain response near the strain at which the twin is activated, where solid lines indicate simulated data, and dashed lines estimate the stress response assuming instantaneous twin activation. Table 15 Change in macroscopic tensile stress during twin activation. ϵ twin (%)

w (mm)

∆σ (MPa)

∆σ (%)

0.5

0.0025 0.005 0.01

−2.0 −4.3 −7.1

−0.31 −0.67 −1.10

1.0

0.0025 0.005 0.01

−8.1 −11.4 −15.1

−0.76 −1.04 −1.41

region – which in this specific case aligns the crystal to a more compliant orientation – and the elastic unloading that occurs upon twin activation, described above. 8. Summary and outlook In this study, an existing crystal plasticity finite element framework is extended to include discrete deformation twinning in full-field three-dimensional simulations by means of pre-discretization of grains into lamellar domains that represent potential twin regions. When deformation twinning is activated, velocities are imposed on the nodes that define a designated lamella such that the lamella is mapped to the expected position due to twinning. Concurrently, the crystal lattice is reoriented to the twin orientation. This deformation is forced to occur over a very brief duration such that the shear rate associated with a twinning event is greater than shear rates commonly witnessed due to crystallographic slip. In short, this framework affords the ability to activate discrete deformation twinning at desired points in deformation, and the mechanical response to this event is calculated. The framework is demonstrated with an example in which a twin is inserted in a chosen grain of a titanium polycrystal. Simulations are performed on polycrystalline aggregates, where a twin is activated in a crystal at the center of the aggregate at two arbitrary points in deformation. Three different twin widths, and two different points of twin activation in the loading history are considered. It is shown that the influence of twin thickness and the point of twin activation can be studied, not only in terms of the local stresses in and around the twin, but also in terms of the work. Results where twinning is activated in the elastic regime indicate that large stress concentrations form at the twin tips in neighboring grains and the volume of the regions effected by the twin motion in neighboring grains increases as the twin size increases. These stress concentrations are not apparent when twinning is activated in the plastic regime, as the stress magnitude is constrained by the single crystal yield surface. The macroscopic

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

29

response is similarly altered. The activation of twinning leads to a decrease in the macroscopic tensile stress due to the reorientation of the crystal lattice within the twin, and the elastic unloading within the polycrystal due to the imposed motion of the twin. The work necessary to activate a twin is calculated to understand the energy associated with a twinning event. For thin twins in the elastic regime, the calculated plastic work is less than the total work. This indicates a negative contribution from the work due to elasticity — that is, some portion of the volume is elastically unloading, indicating that twinning causes some elastic relaxation within the polycrystal. As twin size increases, both the total work and plastic work necessary to activate a twin increase. For thick twins, the total work is found to be greater than the plastic work, indicating net elastic loading throughout the aggregate, rather than relaxation. In the plastic regime, plastic work is greater than the total work for every twin thickness. As described in this study, the framework considers no nucleation criteria or growth criteria. Twins are activated at pre-determined points in the load history, and at pre-determined spatial locations. However, implementation of various nucleation or growth criteria is possible in this framework, and it is envisioned that it will be used for analyses of such criteria. The type of simulation discussed in this study is expected to closely align with high energy X-ray diffraction (HEXD) experiments. As the technology associated with X-ray science progresses, there is the potential to develop methods to measure the development of the strain (stress) field around a twin during in-situ loading. As stated in the introduction, so-called near field HEXD has the potential non-destructively elucidate three-dimensional grain (twin) morphology, while far-field HEXD has been demonstrated to measure the average strain state of grains. One could imagine an experiment where deformation twins are observed in both the near field and far-field measurements, and the average stress states measured. Data concerning twin activation (parent orientation, grain neighborhood, morphology, and strain at activation) could be used to inform a representative simulation, and output data (average stress states) could be compared to narrow in on uncertainty in the prediction of the stress fields. Acknowledgments MK was supported by the Cornell High Energy Synchrotron Source (CHESS), United States of America through this study, under NSF award DMR-1332208. PD received support from the Office of Naval Research, United States of America, under Grant N00014-16-1-2982. Many thanks to Romain Quey for his insightful discussion and helpful suggestions. References [1] G. Lütjering, J. Williams, Titanium, second ed., in: Engineering Materials and Processes, Springer, Berlin; New York, 2007. [2] W. Hosford, The Mechanics of Crystals and Textured Polycrystals, in: The Oxford Engineering Science Series, vol. 32, Oxford University Press, New York, 1993. [3] U.F. Kocks, C.N. Tomé, H. Wenk, Texture and Anisotropy: Preferred Orientations in Polycrystals and their Effect on Materials Properties, first ed., Cambridge Univ. Press, Cambridge, 2000. [4] A. Cottrell, Dislocations and Plastic Flow in Crystals, Oxford University Press, 1953. [5] J. Nye, Physical Properties of Crystals: Their Representation by Tensors and Matrices, Clarendon Press ; Oxford University Press, 1984. [6] R. Ogden, Non-linear Elastic Deformations, Dover Publications, 1997. [7] A. Bower, Applied Mechanics of Solids, CRC Press, 2010. [8] S. Wong, P. Dawson, Influence of directional strength-to-stiffness on the elastic-plastic transition of fcc polycrystals under uniaxial tensile loading, Acta Mater. 58 (5) (2010) 1658–1678, http://dx.doi.org/10.1016/j.actamat.2009.11.009. [9] M. Obstalecki, S. Wong, P. Dawson, M. Miller, Quantitative analysis of crystal scale deformation heterogeneity during cyclic plasticity using high-energy X-ray diffraction and finite-element simulation, Acta Mater. 75 (2014) 259–272, http://dx.doi.org/10.1016/j.actamat. 2014.04.059. [10] S. Wong, M. Obstalecki, M. Miller, P. Dawson, Stress and deformation heterogeneity in individual grains within polycrystals subjected to fully reversed cyclic loading, J. Mech. Phys. Solids 79 (2015) 157–185, http://dx.doi.org/10.1016/j.jmps.2015.03.010. [11] F. Roters, P. Eisenlohr, L. Hantcherli, D. Tjahjanto, T. Bieler, D. Raabe, Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: Theory, experiments, applications, Acta Mater. 58 (4) (2010) 1152–1211, http://dx.doi.org/10.1016/j.actamat.2009.10.058. [12] E. Wielewski, D. Boyce, J.-S. Park, M. Miller, P. Dawson, A methodology to determine the elastic moduli of crystals by matching experimental and simulated lattice strain pole figures using discrete harmonics, Acta Mater. 126 (2017) 469–480, http: //dx.doi.org/10.1016/j.actamat.2016.12.026.

30

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

[13] P. Dawson, D. Boyce, J.-S. Park, E. Wielewski, M. Miller, Determining the strengths of HCP slip systems using harmonic analyses of lattice strain distributions, Acta Mater. 144 (2018) 92–106, http://dx.doi.org/10.1016/j.actamat.2017.10.032. [14] I. Jones, W. Hutchinson, Stress-state dependence of slip in Titanium-6Al-4V and other H.C.P. metals, Acta Mater. 23 (9) (1981) 951–968, http://dx.doi.org/10.1016/0001-6160(81)90049-3. [15] X. Song, S.Y. Zhang, D. Dini, A.M. Korsunsky, Finite element modelling and diffraction measurement of elastic strains during tensile deformation of HCP polycrystals, Comput. Mater. Sci. 44 (1) (2008) 131–137, http://dx.doi.org/10.1016/j.commatsci.2008.01.043. [16] F. Bridier, D. McDowell, P. Villechaise, J. Mendez, Crystal plasticity modeling of slip activity in Ti-6Al-4V under high cycle fatigue loading, Int. J. Plast. 25 (6) (2009) 1457–1485, http://dx.doi.org/10.1016/j.ijplas.2008.08.004. [17] V. Hasija, S. Ghosh, M. Mills, D. Joseph, Deformation and creep modeling in polycrystalline Ti-6Al alloys, Acta Mater. 51 (15) (2003) 4533–4549, http://dx.doi.org/10.1016/s1359-6454(03)00289-1. [18] J.C. Williams, R.G. Baggerly, N.E. Paton, Deformation behavior of HCP Ti-Al alloy single crystals, Metall. Mater. Trans. A 33 (3) (2002) 837–850, http://dx.doi.org/10.1007/s11661-002-0153-y. [19] Y. Yang, L. Wang, C. Zambaldi, P. Eisenlohr, R. Barabash, W. Liu, M.R. Stoudt, M.A. Crimp, T.R. Bieler, Characterization and modeling of heterogeneous deformation in commercial purity titanium, JOM 63 (9) (2011) 66–73, http://dx.doi.org/10.1007/s11837-011-0161-8. [20] J. Christian, S. Mahajan, Deformation twinning, Prog. Mater. Sci. 39 (1–2) (1995) 1–157, http://dx.doi.org/10.1016/0079-6425(94)000077. [21] R. Reed-Hill, Physical Metallurgy Principles, Van Nostrand Reinhold Inc., U.S., 1967. [22] H.J. Frost, M.F. Ashby, Deformation-Mechanism Maps: The Plasticity and Creep of Metals and Ceramics, Pergamon Press, 1982. [23] R.W. Hertzberg, Deformation and Fracture Mechanics of Engineering Materials, Wiley, 1983. [24] G. Dieter, Mechanical Metallurgy, McGraw-Hill Education, 1986. [25] D.C. Stouffer, L.T. Dame, Inelastic Deformation of Metals: Models, Mechanical Properties, and Metallurgy, Wiley-Interscience, 1996. [26] S. Wong, M. Madivala, U. Prahl, F. Roters, D. Raabe, A crystal plasticity model for twinning- and transformation-induced plasticity, Acta Mater. 118 (2016) 140–151, http://dx.doi.org/10.1016/j.actamat.2016.07.032. [27] W.H. Press, Numerical Recipes in C: The Art of Scientific Computing, second ed., Cambridge University Press, 1992. [28] P. Van Houtte, Simulation of the rolling and shear texture of brass by the taylor theory adapted for mechanical twinning, Acta Metall. 26 (4) (1978) 591–604, http://dx.doi.org/10.1016/0001-6160(78)90111-6. [29] C. Tomé, R. Lebensohn, U. Kocks, A model for texture development dominated by deformation twinning: application to zirconium alloys, Acta Metall. Mater. 39 (11) (1991) 2667–2680, http://dx.doi.org/10.1016/0956-7151(91)90083-d. [30] S. Kalidindi, Incorporation of deformation twinning in crystal plasticity models, J. Mech. Phys. Solids 46 (2) (1998) 267–290, http://dx.doi.org/10.1016/s0022-5096(97)00051-3. [31] H. Abdolvand, M. Daymond, C. Mareau, Incorporation of twinning into a crystal plasticity finite element model: Evolution of lattice strains and texture in zircaloy-2, Int. J. Plast. 27 (11) (2011) 1721–1738, http://dx.doi.org/10.1016/j.ijplas.2011.04.005. [32] M. Lindroos, G. Cailletaud, A. Laukkanen, V.-T. Kuokkala, Crystal plasticity modeling and characterization of the deformation twinning and strain hardening in hadfield steels, Mater. Sci. Eng. A 720 (2018) 145–159, http://dx.doi.org/10.1016/j.msea.2018.02.028. [33] R. Quey, P. Dawson, F. Barbe, Large-scale 3D random polycrystals for the finite element method: Generation, meshing and remeshing, Comput. Methods Appl. Mech. Engrg. 200 (17–20) (2011) 1729–1745, http://dx.doi.org/10.1016/j.cma.2011.01.002. [34] R. Quey, L. Renversade, Optimal polyhedral description of 3d polycrystals: Method and application to statistical and synchrotron x-ray diffraction data, Comput. Methods Appl. Mech. Engrg. 330 (2018) 308–333, http://dx.doi.org/10.1016/j.cma.2017.10.029. [35] M. Kasemer, R. Quey, P. Dawson, The influence of mechanical constraints introduced by β annealed microstructures on the yield strength and ductility of Ti-6Al-4V, J. Mech. Phys. Solids 103C (2017) 179–198, http://dx.doi.org/10.1016/j.jmps.2017.03.013. [36] R. Quey, Neper: polycrystal generation and meshing (version 3.0), http://neper.sourceforge.net (2016). [37] E. Marin, P. Dawson, On modelling the elasto-viscoplastic response of metals using polycrystal plasticity, Comput. Methods Appl. Mech. Engrg. 165 (1–4) (1998) 1–21, http://dx.doi.org/10.1016/S0045-7825(98)00034-6. [38] E. Marin, P. Dawson, Elastoplastic finite element analyses of metal deformations using polycrystal constitutive models, Comput. Methods Appl. Mech. Engrg. 165 (1–4) (1998) 23–41, http://dx.doi.org/10.1016/S0045-7825(98)00033-4. [39] M. Kasemer, M. Echlin, J. Stinville, T. Pollock, P. Dawson, On slip initiation in equiaxed α/β Ti-6Al-4V, Acta Mater. 136 (2017) 288–302, http://dx.doi.org/10.1016/j.actamat.2017.06.059. [40] R. Asaro, A. Needleman, Texture development and strain hardening in rate dependent polycrystals, Acta Metall. 33 (6) (1985) 923–953, http://dx.doi.org/10.1016/0001-6160(85)90188-9. [41] N. Barton, P. Dawson, Lattice misorientations in titanium alloys: modeling the origins of defects, Int. J. Forming Process. 5 (2-3-4) (2002) 189–201. [42] C. Tomé, R. Lebensohn, U. Kocks, A model for texture development dominated by deformation twinning: application to zirconium alloys, Acta Metall. Mater. 39 (11) (1991) 2667–2680, http://dx.doi.org/10.1016/0956-7151(91)90083-d. [43] R. Lebensohn, C. Tomé, A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals: Application to zirconium alloys, Acta Metall. Mater. 41 (9) (1993) 2611–2624, http://dx.doi.org/10.1016/09567151(93)90130-k. [44] M.A. Kumar, I. Beyerlein, C. Tomé, A measure of plastic anisotropy for hexagonal close packed metals: Application to alloying effects on the formability of mg, J. Alloys Compd. 695 (2017) 1488–1497, http://dx.doi.org/10.1016/j.jallcom.2016.10.287. [45] T. Heo, Y. Wang, S. Bhattacharya, X. Sun, S. Hu, L. Chen, A phase-field model for deformation twinning, Phil. Mag. Lett. 91 (2) (2011) 110–121, http://dx.doi.org/10.1080/09500839.2010.537284. [46] M.A. Kumar, I. Beyerlein, R. Lebensohn, C. Tomé, Role of alloying elements on twin growth and twin transmission in magnesium alloys, Mater. Sci. Eng. A 706 (2017) 295–303, http://dx.doi.org/10.1016/j.msea.2017.08.084.

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

31

[47] C. Liu, P. Shanthraj, M. Diehl, F. Roters, S. Dong, J. Dong, W. Ding, D. Raabe, An integrated crystal plasticity–phase field model for spatially resolved twin nucleation, propagation, and growth in hexagonal materials, Int. J. Plast. 106 (2018) 203–227, http://dx.doi.org/10.1016/j.ijplas.2018.03.009. [48] M. Ardeljan, R. McCabe, I. Beyerlein, M. Knezevic, Explicit incorporation of deformation twins into crystal plasticity finite element models, Comput. Methods Appl. Mech. Engrg. 295 (2015) 396–413, http://dx.doi.org/10.1016/j.cma.2015.07.003. [49] M. Ardeljan, I.J. Beyerlein, M. Knezevic, Effect of dislocation density-twin interactions on twin growth in AZ31 as revealed by explicit crystal plasticity finite element modeling, Int. J. Plast. 99 (2017) 81–101, http://dx.doi.org/10.1016/j.ijplas.2017.09.002. [50] M. Ardeljan, M. Knezevic, Explicit modeling of double twinning in AZ31 using crystal plasticity finite elements for predicting the mechanical fields for twin variant selection and fracture analyses, Acta Mater. 157 (2018) 339–354, http://dx.doi.org/10.1016/j.actamat. 2018.07.045. [51] J. Cheng, S. Ghosh, Crystal plasticity finite element modeling of discrete twin evolution in polycrystalline magnesium, J. Mech. Phys. Solids 99 (2017) 512–538, http://dx.doi.org/10.1016/j.jmps.2016.12.008. [52] J. Cheng, J. Shen, R. Mishra, S. Ghosh, Discrete twin evolution in mg alloys using a novel crystal plasticity finite element model, Acta Mater. 149 (2018) 142–153, http://dx.doi.org/10.1016/j.actamat.2018.02.032. [53] S. Meftah, F. Barbe, L. Taleb, F. Sidoroff, Parametric numerical simulations of TRIP and its interaction with classical plasticity in martensitic transformation, Eur. J. Mech. A Solids 26 (4) (2007) 688–700, http://dx.doi.org/10.1016/j.euromechsol.2006.10.004. [54] T. Ungár, M. Glavicic, L. Balogh, K. Nyilas, A. Salem, G. Ribárik, S. Semiatin, The use of x-ray diffraction to determine slip and twinning activity in commercial-purity (CP) titanium, Mater. Sci. Eng. A 493 (1–2) (2008) 79–85, http://dx.doi.org/10.1016/j.msea. 2007.06.096. [55] L. Wang, J. Lind, H. Phukan, P. Kenesei, J.-S. Park, R. Suter, A. Beaudoin, T. Bieler, Mechanical twinning and detwinning in pure ti during loading and unloading - an in situ high-energy x-ray diffraction microscopy study, Scr. Mater. 92 (2014) 35–38, http://dx.doi.org/10.1016/j.scriptamat.2014.08.008. [56] T. Ungár, Microstructural parameters from x-ray diffraction peak broadening, Scr. Mater. 51 (8) (2004) 777–781, http://dx.doi.org/10. 1016/j.scriptamat.2004.05.007. [57] W. Ludwig, S. Schmidt, E.M. Lauridsen, H.F. Poulsen, X-ray diffraction contrast tomography: a novel technique for threedimensional grain mapping of polycrystals. i. direct beam case, J. Appl. Crystallogr. 41 (2) (2008) 302–309, http://dx.doi.org/10. 1107/s0021889808001684. [58] U. Lienert, S.F. Li, C.M. Hefferan, J. Lind, R.M. Suter, J.V. Bernier, N.R. Barton, M.C. Brandes, M.J. Mills, M.P. Miller, B. Jakobsen, W. Pantleon, High-energy diffraction microscopy at the advanced photon source, JOM 63 (7) (2011) 70–77, http: //dx.doi.org/10.1007/s11837-011-0116-0. [59] H. Abdolvand, M. Majkut, J. Oddershede, J. Wright, M. Daymond, Study of 3-d stress development in parent and twin pairs of a hexagonal close-packed polycrystal: Part II – crystal plasticity finite element modeling, Acta Mater. 93 (2015) 235–245, http://dx.doi.org/10.1016/j.actamat.2015.04.025. [60] E. Schmid, W. Boas, Kristallplastizität mit besonderer Berücksichtigung der Metalle, Springer-Verlag, Berlin, 1935. [61] S. Myagchilov, P. Dawson, Evolution of texture in aggregates of crystals exhibiting both slip and twinning, Modelling Simulation Mater. Sci. Eng. 7 (6) (1999) 975–1004, http://dx.doi.org/10.1088/0965-0393/7/6/305. [62] P.R. Dawson, D.E. Boyce, FEpX Finite Element Polycrystals: Theory, Finite Element Formulation, Numerical Implementation and Illustrative Examples, arXiv e-prints (2015) arXiv:1504.03296arXiv:1504.03296. [63] A. Ghaderi, M. Barnett, Sensitivity of deformation twinning to grain size in titanium and magnesium, Acta Mater. 59 (20) (2011) 7824–7839, http://dx.doi.org/10.1016/j.actamat.2011.09.018. [64] F. Frank, Orientation mapping, MRS Bull. 13 (03) (1988) 24–31, http://dx.doi.org/10.1557/S0883769400066112. [65] Z. Zhang, T.-S. Jun, T.B. Britton, F.P. Dunne, Intrinsic anisotropy of strain rate sensitivity in single crystal alpha titanium, Acta Mater. 118 (2016) 317–330, http://dx.doi.org/10.1016/j.actamat.2016.07.044. [66] Z. Zheng, S. Waheed, D.S. Balint, F.P. Dunne, Slip transfer across phase boundaries in dual phase titanium alloys and the effect on strain rate sensitivity, Int. J. Plast. 104 (2018) 23–38, http://dx.doi.org/10.1016/j.ijplas.2018.01.011. [67] A. Salem, S. Kalidindi, S. Semiatin, Strain hardening due to deformation twinning in α-titanium: Constitutive relations and crystal-plasticity modeling, Acta Mater. 53 (12) (2005) 3495–3502, http://dx.doi.org/10.1016/j.actamat.2005.04.014. [68] M. Niewczas, Lattice correspondence during twinning in hexagonal close-packed crystals, Acta Mater. 58 (17) (2010) 5848–5857, http://dx.doi.org/10.1016/j.actamat.2010.06.059. [69] A. Akhtar, Basal slip and twinning in α-titanium single crystals, Metall. Trans. A 6 (5) (1975) 1105–1113, http://dx.doi.org/10.1007/ bf02661366. [70] D. Rowenhorst, A.D. Rollett, G.S. Rohrer, M. Groeber, M. Jackson, P.J. Konijnenberg, M.D. Graef, Consistent representations of and conversions between 3d rotations, Modelling Simulation Mater. Sci. Eng. 23 (8) (2015) 083501, http://dx.doi.org/10.1088/09650393/23/8/083501. [71] H. Abdolvand, A. Wilkinson, Assessment of residual stress fields at deformation twin tips and the surrounding environments, Acta Mater. 105 (2016) 219–231, http://dx.doi.org/10.1016/j.actamat.2015.11.036. [72] Z. Zhang, D. Lunt, H. Abdolvand, A. Wilkinson, M. Preuss, F. Dunne, Quantitative investigation of micro slip and localization in polycrystalline materials under uniaxial tension, Int. J. Plast. 108 (2018) 88–106, http://dx.doi.org/10.1016/j.ijplas.2018.04.014. [73] K. Kapoor, Y. Yoo, T. Book, J. Kacher, M. Sangid, Incorporating grain-level residual stresses and validating a crystal plasticity model of a two-phase ti-6al-4v alloy produced via additive manufacturing, J. Mech. Phys. Solids 121 (2018) 447–462, http: //dx.doi.org/10.1016/j.jmps.2018.07.025. [74] Z. Zhang, F. Dunne, Phase morphology, variants and crystallography of alloy microstructures in cold dwell fatigue, Int. J. Fatigue 113 (2018) 324–334, http://dx.doi.org/10.1016/j.ijfatigue.2018.03.030.

32

M. Kasemer and P. Dawson / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112653

[75] K. Chatterjee, M. Echlin, M. Kasemer, P. Callahan, T. Pollock, P. Dawson, Prediction of tensile stiffness and strength of Ti-6Al-4V using instantiated volume elements and crystal plasticity, Acta Mater. 157 (2018) 21–32, http://dx.doi.org/10.1016/j.actamat.2018.07.011. [76] U. Ayachit, The ParaView guide: Updated for ParaView Version 4.3, Kitware, Clifton Park, New York, 2015. [77] M.A. Kumar, I.J. Beyerlein, R.A. Lebensohn, C.N. Tomé, Modeling the effect of neighboring grains on twin growth in HCP polycrystals, Modelling Simulation Mater. Sci. Eng. 25 (6) (2017) 064007, http://dx.doi.org/10.1088/1361-651x/aa7bbb. [78] I. Chelladurai, D. Adams, D.T. Fullwood, M.P. Miles, S. Niezgoda, I.J. Beyerlein, M. Knezevic, Modeling of trans-grain twin transmission in AZ31 via a neighborhood-based viscoplastic self-consistent model, Int. J. Plast. 117 (2019) 21–32, http://dx.doi.org/ 10.1016/j.ijplas.2018.03.012.