J Mater Sci I N T E R F AC E Behavior B E H A V I OR Interface
Mesoscale models of interface mechanics in crystalline solids: a review J. D. Clayton1,2,3,* 1
Department of Civil Engineering and Engineering Mechanics, Columbia University, New York, NY 10027, USA A. James Clark School of Engineering, University of Maryland, College Park, MD 20742, USA 3 Impact Physics, US ARL, Aberdeen, MD 21005-5066, USA 2
Received: 14 June 2017
ABSTRACT
Accepted: 18 September 2017
Theoretical and computational methods for representing mechanical behaviors of crystalline materials in the vicinity of planar interfaces are examined and compared. Emphasis is on continuum-type resolutions of microstructures at the nanometer and micrometer levels, i.e., mesoscale models. Grain boundary interfaces are considered first, with classes of models encompassing sharp interface, continuum defect (i.e., dislocation and disclination), and diffuse interface types. Twin boundaries are reviewed next, considering sharp interface and diffuse interface (e.g., phase field) models as well as pseudo-slip crystal plasticity approaches to deformation twinning. Several classes of models for evolving failure interfaces, i.e., fracture surfaces, in single crystals and polycrystals are then critically summarized, including cohesive zone approaches, continuum damage theories, and diffuse interface models. Important characteristics of compared classes of models for a given physical behavior include complexity, generality/flexibility, and predictive capability versus number of free or calibrated parameters.
Ó
Springer Science+Business
Media, LLC 2017
Introduction Interfaces may have a profound effect on the physical responses of solid materials subjected to various external stimuli or environmental conditions. This article is focused on crystalline solid materials subjected to mechanical loading [1–3]. In particular, material classes considered herein may be single
crystals or polycrystals, comprised of grains of metallic, ceramic, molecular, and/or mineralogical origin. In order to maintain a reasonable scope, not considered in this article are polymers and polymer composites, though it is acknowledged that interfacial mechanics is often crucial to the local and global responses of these classes of materials as well [4, 5]. An interface is defined herein as a planar boundary between two regions of a solid across which physical
Invited review article for Special Issue of Journal of Materials Science.
Address correspondence to E-mail:
[email protected];
[email protected]
DOI 10.1007/s10853-017-1596-2
J Mater Sci
properties and/or response characteristics vary substantially. Local behaviors of the material in the immediate vicinity of such an interface may also vary from that of bulk regions of crystal far from the boundary on either side. Examples of interface types that may exist prior to deformation or mechanical loading include grain [6] and subgrain boundaries, twin boundaries, and phase boundaries. Those that are often induced by loading include fracture surfaces [7], stacking faults, and shear bands [8]. Grain and twin boundaries may also be deformation induced, as in cases of grain subdivision during plastic deformation [9] or deformation twinning [10, 11]. This paper is focused on modeling techniques, including theory and methods of solution, both analytical and numerical. Experiments are mentioned at times in supporting discussion, e.g., in the context of model validation, but this paper does not critically examine experimental methods. The scope is limited to continuum models that resolve microstructural features, labeled here as mesoscale models. Length scales resolved span approximately tens of nanometers to millimeters, with micrometers the primary scale of most representations discussed herein. Not addressed are atomic scale theories [12, 13], or coupled atomic-continuum models [14–17] which tend to focus on somewhat finer scales of resolution. General classes of models include sharp interface theories, whereby jumps in properties and/or field variables exist across internal boundaries, and diffuse interface theories, wherein some smoothing or continuous interpolation of properties and/or response functions is invoked across an interface [18]. The overall scope of this paper is thought to be unique in terms of the ensemble of topics and their comparisons. An attempt to sufficiently mention relevant works in each area is made, but inclusion of even a modest fraction of the vast number of historical and recent articles on any broad subject is currently impossible for a review paper of any reasonably finite length, especially when one considers prolific publishing practices in contemporary fundamental sciences. Later in this work, for each kind of represented physics—grain boundaries, twin boundaries, failure boundaries—the corresponding classes of models are compared with regard to flexibility, complexity, and predictive capability. Flexibility, i.e., generality, refers to the ability of a model to represent a breadth of physical behaviors without adjustment of its fundamental governing equations. Complexity of a theory
is rather self-explanatory, with availability of analytical solutions and ease of numerical implementation both tending to decrease with increasing complexity. The viewpoint adopted here is that a given theory or model is considered more predictive than another if it equally or better represents observed physics with fewer calibrated parameters [19–21]. This paper is organized as follows, noting that as is the case with nearly any review, the present subject matter deals with topics of the author’s own experience. The ‘‘Grain boundaries’’ section addresses grain boundaries in polycrystals, including sharp interface, continuum defect, and diffuse interface or phase field representations. The ‘‘Twinning’’ section covers twin boundaries, again considering sharp interface and phase field models as well as pseudo-slip theories for deformation twinning. The ‘‘Failure and localization’’ section describes fracture and localization models in single crystals and polycrystals: cohesive zone methods, continuum damage models, and regularized failure models such as those of phase field type. Fractures may be transgranular or intergranular, with spall fractures potentially falling into either category. The ‘‘Discussion’’ section summarizes major features of the approaches including their strengths and weaknesses according to the author’s viewpoint. Although this review is biased more toward general themes than intensive derivations, inclusion of key mathematical relations for theories and methods of solution is necessary for the intended rigorous evaluation. Notation of continuum physics is invoked in such instances, where vectors and tensors are written in bold italic font, and scalars and scalar components in italic font. Occasionally, the index notation is used for clarity, with summation implied over repeated indices. Further details regarding notation will be clear from context.
Grain boundaries Prior to discussion of specific classes of models in the context of mechanical response, some universal features and requisite notation for a grain boundary (GB) are introduced. Useful references for further study are [22, 23]. Superscripts in parentheses are used to denote grain numbers in a polycrystal, where such a number runs from 1 to the total number of crystals within the
J Mater Sci
body. Denote by XðgÞ the volume occupied by grain g, prior to deformation. Let X denote the position vector of a material point in the body, again prior to deformation induced by loading, and referred to a fixed Cartesian reference frame. Let oXðgÞ denote the external boundary of a crystal, which may include free surfaces and/or grain boundaries (GBs). The boundary between two grains with numbers g1 and g2 is denoted by oXðg1;g2Þ , and so forth. Bravais lattice vectors ai , where i ¼ 1; 2; 3, are assigned to each single crystal and are assumed uniform in each crystal volume, i.e., ðgÞ
ai ðXÞ ¼ ai
8X 2 XðgÞ :
ð2:1Þ
By definition, the Bravais lattice is discontinuous across (misoriented) grain boundaries. Let a0i denote a particular set of the Bravais lattice vectors referred to a given coordinate frame. Then, the orientation of the crystal lattice in any grain is related to these vectors via the orientation matrix GðgÞ : ðgÞ
ai
¼ GðgÞ a0i :
ð2:2Þ
The misorientation matrix for any two crystals with grain numbers g1 and g2 is defined as M ðg1;g2Þ ¼ ðGðg1Þ Þ1 Gðg2Þ :
ð2:3Þ
Orientation matrices and misorientation matrices are all proper orthogonal, i.e., each has equal transpose and inverse, with a determinant of value þ1. Such matrices can be represented in terms of an angle u and axis of rotation r of unit magnitude. For a generic misorientation matrix M with entries MIJ (I; J ¼ 1; 2; 3), the angle and components of the axis are calculated as follows [23]: 1 cos u ¼ ðM11 þ M22 þ M33 1Þ; 2 r1 ¼ M23 M32 ;
r2 ¼ M31 M13 ;
ð2:4Þ r3 ¼ M12 M21 : ð2:5Þ
Often later in this article, the term misorientation will be used to refer to the magnitude juj, where it is understood that a complete description of the geometric relation between orientations of the two lattices of the same crystal structure requires the axis r as well. Furthermore, it is noted that a complete description of a planar interface requires five independent scalar parameters; the misorientation matrix or angle/axis pair each supplies only three. For
example, a complete characterization of a GB orientation can be achieved via specification of the direction cosines of the boundary plane in the coordinate systems of each neighboring grain along with the angle of twist of both plane stacks normal to the boundary plane [23]. The coincident site lattice (CSL) model of GBs [6, 24] will be used in some later descriptions. For a certain misorientation across a GB, or between two interpenetrating lattices, some fraction of lattice sites will coincide, forming a periodic sublattice. The CSL parameter R is defined as the reciprocal density of coinciding sites. Boundaries associated with low values of R are of high interest since special physical properties are often a result [23]. The use of the grain boundary character distribution to design polycrystalline materials with optimum physical properties has been suggested [25]. However, a low R value does not always correlate with GB strength [26]. The maximum deviation from an exact CSL usually is taken to correspond to the angular limit for a lowangle boundary: juj ¼ 15 . Thus, a low-angle boundary is described as R1. Annealing twin boundaries are a subset of R3 boundaries. An empirical relationship between maximum angular deviation um of an arbitrary GB from any exact CSL with value R is um ¼ u0 R1=2 , where u0 is the aforementioned maximum deviation, typically 15 [6, 23]. The jump in any field variable or material property AðXÞ, which could be a scalar, vector, or higher-order tensor, is defined as the difference between its values at two corresponding locations X ð1Þ and X ð2Þ : sAðXð1Þ ; Xð2Þ Þt ¼ AðX ð2Þ Þ AðX ð1Þ Þ:
ð2:6Þ
The jump across a (sharp) grain boundary corresponds to the difference in limiting values of A as X approaches the shared point on boundary from either side. In this work, attention is restricted to coherent interfaces in the undeformed reference configuration, meaning sXt ¼ 0. Later, in the ‘‘Failure and localization’’ section, surfaces of discontinuity such as fractured interfaces will be discussed, but these are assumed to be induced by loading or deformation from a coherent reference state.
Sharp interface models In what is termed herein as a sharp interface model of a GB, properties of crystals in the immediate vicinity of the boundary are identical to those in the interior
J Mater Sci
of each crystal, far from the interface. Therefore, according to (2.3) and (2.6), the lattice vectors will demonstrate a discrete jump across the GB for any nontrivial misorientation matrix between the two grains: ðg2Þ
sai t ¼ ai
ðg1Þ
ai
¼ ½Gðg2Þ Gðg1Þ a0i
¼ Gðg1Þ ½M ðg1;g2Þ 1a0i
on oXðg1;g2Þ :
ð2:7Þ
The identity tensor is denoted by 1. The difference in orientation of the lattice across the boundary can induce various local physical responses when the aggregate is subjected to far-field mechanical loading, temperature change, and so forth. A geometric representation of a sharp interface discretization of a polycrystal is shown on the left side of Fig. 1, where individual crystals are polyhedral shaped and GBs are faceted (flat) planar surfaces. Such representations are characteristic of modern finite element simulations of crystal elasticity [27, 28] and crystal plasticity [26, 29]. This work will deal with potentially large strains and rotations, as may occur during deformation of ductile metals [3, 26], and even in ceramics and minerals when loading is predominantly compressive [31, 32]. Governing equations from finite elasticity and plasticity of crystals are now reviewed to lend context to the discussion of interfacial mechanics. Only essential theoretical relations are provided, and Cartesian coordinates are implied when index notation is invoked. For a more comprehensive treatment of finite crystal elastoplasticity that also encompasses curvilinear coordinate representations,
see [3]. Two other useful references, albeit primarily limited to finite anisotropic elasticity, are [33, 34]. Spatial coordinates of the deformed solid are related to material coordinates by the time-dependent motion x ¼ xðX; tÞ ¼ X þ uðX; tÞ;
ð2:8Þ
with u the displacement vector. The deformation gradient is, with r0 ðÞ the referential gradient, F ¼ r0 x ¼ FE FP :
ð2:9Þ
Thermoelastic and plastic deformations are FE and FP , respectively. The plastic velocity gradient is the sum of contributions of slip rates c_a , where the superscript denotes a slip system for dislocation glide with direction sa and plane normal ma : X LP ¼ F_ P FP1 ¼ c_a sa ma : ð2:10Þ a
Here the slip direction and slip plane normal are orthogonal and of unit length, i.e., are those of the crystal lattice prior to thermoelastic deformation. The thermoelastic strain used in standard crystal hyperelasticity [33] is the Green strain tensor: i 1h 1 E ¼ ðFE ÞT FE 1 ; EIJ ¼ ðFEkI FEkJ dIJ Þ: 2 2 ð2:11Þ The thermoelastic volume change is measured by J E ¼ det FE . Cauchy stress r is related to elastic second Piola–Kirchhoff stress S via
Figure 1 Geometric rendering of polycrystal with grain boundaries (GBs) represented as sharp interfaces: without secondary GB phase (left); with secondary GB phase (right) [30].
J Mater Sci
r¼
1 E F SðFE ÞT ; JE
rij ¼
1 E F SKL FEjL : J E iK
ð2:12Þ
CðXÞ ¼ C½ai ðXÞ ¼ C½GðXÞ; sa ðXÞ ¼ sa ½ai ðXÞ ¼ sa ½GðXÞ; ma ðXÞ ¼ ma ½ai ðXÞ ¼ ma ½GðXÞ: ð2:17Þ
The thermoelastic stress–strain relation is SIJ ¼ CIJKL EKL þ þ
1 CIJKLMN EKL EMN 2!
1 CIJKLMNPQ EKL EMN EPQ þ vIJ DT ; 3! ð2:13Þ
where CIJKL are isothermal elastic constants of second and higher orders, DT is temperature change measured from a reference state, and vIJ are thermal stress coefficients related to elastic constants and thermal expansion coefficients via vIJ ¼ CIJKL aKL . A free energy function per unit volume in the thermoelastically unloaded configuration is w ¼ wðE; T; fngÞ;
ð2:14Þ
with fng a set of internal state variables that affect the energy stored in the crystal, e.g., dislocation density. The viscoplastic flow rule for slip rates is of the general form c_a ¼ c_a ðsa ; T; fngÞ;
sa ¼ J E r : ½FE sa ðFE ÞT ma : ð2:15Þ
The resolved shear stress acting on a system is sa . The local balance of energy, in the absence of point heat sources, can be cast as the following temperature rate equation: X ow o2 w _ Tv : E_ cT_ ¼ T sa c_a fng ofng oTofng a X ¼b sa c_a Tv : E_ þ r ðK rTÞ; þ r ðK rTÞ a
ð2:16Þ where c is the specific heat at constant thermoelastic strain, the rightmost term accounts for heat conduction (see e.g., [3, 21] for details), and b is the fraction of plastic work converted to heat energy, i.e., the fraction of stored energy of cold work is 1 b. Let C denote any of the tensors of elastic moduli in (2.13). These elastic constants, which are the usual second-order type as well as third- and fourth-order constants important in high-pressure applications [31, 33, 35, 36], depend on the orientation of the crystal lattice in the reference configuration, as do the slip directors and slip plane normal vectors:
Similarly, thermal stress coefficients v and thermal conductivity K also depend on crystal orientation in the unloaded state for Lagrangian thermoelasticity. Functional forms of material coefficients for various crystal classes are available in [3, 33, 34]. In the present context of fully coherent GBs, and in the absence of any transgranular separation modes, continuity of displacement and traction in grain interiors and/or along GBs is stated mathematically as suðX; tÞt ¼ 0 , sxðX; tÞt ¼ 0; ð8X 2 XðgÞ and stt ¼ srtn ¼ 0;
8X on oXðg1;g2Þ Þ;
ð8Xðx; tÞon oXðg1;g2Þ Þ:
ð2:18Þ ð2:19Þ
In the traction continuity equation (2.19), n is normal to the surface, and it is assumed that shock waves (i.e., stress jumps in solid dynamics) are absent. Consider a general boundary value problem for a polycrystal, where traction and/or displacement are imposed along the far-field (external) boundary. The continuity equations in (2.18) and (2.19) impose restrictions on the solution strain and stress fields E and r. In the present sharp interface representation, the jump conditions on the lattice vectors in (2.7) lead to jumps in elastic constants and slip director/normal vectors in (2.17). It is these potentially drastic changes in local properties, taken in combination with the continuity constraints for coherent sharp GBs, that may give rise to concentrated stress, strain, and/or slip activity in the vicinity of GBs in solutions to general polycrystal boundary value problems [37]. The propensity for statistically more pronounced stress concentrations at GBs and triple point junctions (i.e., where three GBs intersect) has been predicted in finite deformation crystal plasticity simulations [29, 38]. Lattice mismatch does not necessarily preclude development of continuous bands, across GBs, of homogeneous or localized field variables such as plastic strain, temperature, and/or dislocation density (when modeled as an internal variable depending on cumulative slip), as demonstrated in early [39, 40] and more recent [41] finite element simulations of shear banding in polycrystals. Preferential lattice orientations may promote localized slip bands
J Mater Sci
Figure 2 Predicted strain localization in Al polycrystal at applied tensile strain of 8%: plastic strain, temperature, heat dissipation fraction b, and dislocation density [41].
in isothermal simulations [39, 40], whereas the combination of temperature rise from (2.16) in conjunction with increased slip activity due to thermal activation if included in (2.15) will promote adiabatic shear banding [8]. As is evident from Fig. 2, correlation between dissipation fraction b, dislocation density n, and local strain and temperature concentrations has been predicted in mesoscale crystal plasticity simulations [41]. The earliest finite element simulations of finite (poly)crystal plasticity appeared in the early to mid1980s [42]. Thereafter, transmission and/or blockage of slip at special (CSL R9 and R17b) high-angle boundaries in copper have been studied via computational crystal plasticity and sharp interface representations of GBs [43]; reported results are in general qualitative agreement with experimental trends relating dislocation pile ups and suggested pending fracture mechanisms. Although typically less prone than most metals to plastic deformation from dislocation glide, ceramics and minerals also demonstrate concentrated thermoelastic strain and stress concentrations due to
Figure 3 Tilt boundary (left) comprised of edge dislocations (center) or partial disclination dipoles (right) [3, 53].
lattice mismatch at GBs. A relatively early numerical study of stress concentrations in polycrystals with sharp interface GB representations is [44], which considers elastic anisotropy and anisotropic thermal expansion coefficients and concludes that anisotropy
J Mater Sci
of either property may affect statistical distributions of tractions at GBs. This study will also be referred to later in the ‘‘Failure and Localization’’ section in the context of failure models, since it includes results incorporating cohesive fracture elements at GBs for some simulations. An important physical feature of many industrial ceramics is the presence of secondary, often glassy, phases of impurities at GBs [30, 45]. In a sharp interface model, such phases can be introduced via a thin layer of secondary material between crystals as shown on the right side of Fig. 1. In summary, sharp interface representations of GBs allow for explicit modeling of strain and stress concentrations of field variables near mismatched interfaces. No additional constitutive model parameters are required beyond those invoked for the response of the bulk material in grain interiors. The approach is considered physically realistic for modeling continuum thermoelasticity since thermoelastic properties and hardness tend to vary little from points near the GB to points far in the grain interior [46]. In contrast, properties associated with dissipative phenomena, including dislocation nucleation and motion in ductile crystals, may vary strongly in the vicinity of GBs, as will be discussed more in the ‘‘Continuum dislocation and disclination models’’ section. Furthermore, since no intrinsic length scale exists in classical crystal elasticity and plasticity theories, the sharp interface models that invoke classical bulk constitutive laws are unable to predict size effects such as Hall–Petch-type behavior (i.e., strength increases as grain size decreases). The preceding treatment is focused on stationary grain boundaries. Depending on thermomechanical loading conditions, the mobility of GBs becomes important, for example under conditions pertinent to grain growth, recrystallization, or some kinds of phase transformations. Sharp interface models of surfaces of strain discontinuity and phase boundaries—including geometric, kinematic, thermodynamic, and kinetic aspects—were advanced in the late twentieth century in several seminal works [47–49].
history. As shown on the left side of Fig. 3, a tilt boundary with misorientation magnitude u can be represented in terms of a density of edge dislocations of like sign and Burgers vector magnitude b, with spacing h given by [50] h ¼ b=u;
ðu . u0 ¼ p=12Þ:
ð2:20Þ
To a certain extent, the larger the misorientation, the more closely spaced the dislocations, leading to an overall higher density (line number per unit area). However, the maximum misorientation is limited to around 15 in this model since high-angle grain boundaries are not well represented physically by large dislocation densities. The low-angle boundaries represented in such cases may be present in the undeformed (poly)crystal, or may be induced by severe plastic deformation [9]. Yield behavior has been observed in indentation experiments to vary with distance from GBs in metals, with such dependency varying with orientation mismatch [46]. Boundaries may act as sources or sinks for dislocation generation and interactions. Furthermore, slip transmission or blockage may ensue depending on alignment of slip systems across a GB [43]. Early continuum models accounting for different plastic properties near GBs versus grain interiors treated each crystal as a composite material with a boundary layer of properties differing from those of the bulk [37, 51]. More recent computational models have incorporated rules for GB slip interactions via specialized interfacial finite elements along grain or subgrain boundaries [26, 52]. The GB dislocations concentrated at interfaces between severely deformed crystals, which are required to maintain compatibility of the total deformation gradient field, account for the incompatibility (i.e., failure of integrability or anholonomic character [54–57]) of the thermoelastic and plastic deformations individually in (2.9). The spatial density tensor of geometrically necessary dislocations (GNDs) is, in coordinate free tensor notation [3, 57] : e ¼ T : e ¼ fFE ½rðFE 1 Þg : e; qG ¼ C
ð2:21Þ
Continuum dislocation and disclination models
is the integrable where e is the permutation tensor, C crystal connection whose curvature tensor vanishes identically, and T is the torsion tensor of this con-
Grain boundaries, particularly those of low angle type (i.e., R1), often contain high densities of dislocations depending on the material and processing
nection. The spatial gradient operator is rðÞ. Since the GND tensor in (2.21) does not account for distributions of dislocations whose net Burgers vector vanishes, a scalar density of statistically stored
J Mater Sci
dislocations (SSDs), qS , is also introduced. Let qG be the frame invariant GND tensor mapped to the elastically unloaded intermediate configuration [3, 57–59], a tensor which can equivalently be constructed from the transformation of (2.21) or the material gradient of FP . Then, the state variable list for the crystal plasticity constitutive model entering the free energy and/or flow rule in (2.14) and (2.15) is, for the present class of models, fng ¼ fqG ; qS g:
ð2:22Þ
The flow rule is supplemented with a kinetic equation for the rate of change of qS , which generally tends to increase with cumulative plastic slip. Model parameters are needed to scale the magnitude of the contribution of GNDs to stored energy and hardening behaviors. Similarly, parameters are required to quantify contributions of SSDs to stored energy and slip resistance. A common practice involves use of an effective shear modulus and scalar Burgers vector to provide the correct units for normalization. It is understood that (2.22) must be extended if slip systems harden unequally (i.e., different self- and latent hardening), in which case dislocation densities relevant to each system may become distinct state variables [52]. Inclusion of the GND density tensor in the constitutive model provides a regularizing effect on numerical solutions due to the presence of the gradient term [i.e., rðFE 1 Þ] in (2.21). Size effects can be predicted using this class of models if properly calibrated; for example, such models may depict increased strength under indentation with a decrease in indenter size, or Hall–Petch-type behavior (essentially, smaller ¼ stronger) [26]. In simulations of polycrystals, GNDs tend to build up in the vicinity of grain boundaries, thereby quantifying intergranular incompatibility and physically representing pile ups and increased hardness in such regions [26, 52]. Boundary conditions supplementing those of classical crystal plasticity are required for solution of boundary value problems, e.g., the thermoelastic or plastic deformation gradient may need to be specified over some part of the external boundary of the domain. A recent novel scheme for numerical representation of grain boundaries in polycrystals via smoothed jumps of local lattice rotation in the reference state is described in [60].
A more sophisticated class of elastic–plastic models intended to capture physics of initial and/or evolving boundaries of lattice misorientation in crystals supplements the dislocation (GND and SSD) description with disclinations. Dislocations are fundamental translational defects, while disclinations are fundamental rotational defects, so the latter perhaps more naturally should be used to represent rotations of the lattice across boundaries of misorientation. As shown on the right side of Fig. 3, the mean spacing l of defects in a tilt boundary represented by a density of parallel partial wedge disclination dipoles of strength x and dipole radius r is l ¼ 2rx=u:
ð2:23Þ
The magnitude of misorientation is u, not necessarily restricted to small angles. This disclination description credited to Li [61] is an alternative to the Read– Shockley model of (2.20). With r and x fundamental constant properties of the crystal lattice associated with a particular material (i.e., rotational analogs of b), (2.23) suggests a generally increasing density of disclinations per unit area with increasing misorientation in the vicinity of grain boundaries. In the context of finite deformation continuum field theories of elastic–plastic crystal mechanics, a spatial tensor of disclination density h is defined as [3, 53, 62] 1 h ¼ ðR : eÞ : e; 4
R ¼ RðC; rCÞ;
þ Q: C¼C ð2:24Þ
In differential-geometric language [56, 57, 63], R is the fourth-order curvature tensor that depends on the connection coefficients of C and spatial gradients of these coefficients. The latter total connection C consists of the sum of the crystal connection defined in (2.21) and the micro-rotation variable Q, which physically represents lattice rotation gradients [62]. When Q vanishes, the curvature tensor and the disclination density tensor also vanish. The spatial GND density is constructed as in (2.21), but with C replaced with C. The list of state variables entering the free energy function and flow rule is extended to include disclination density mapped to the intermediate configuration h, i.e., fng ¼ fqG ; qS ; hg:
ð2:25Þ
The list may be further augmented to include a scalar density of statistically stored disclinations with
J Mater Sci
vanishing average Frank vector [3, 53], but such a description is not essential in the present context. Continuum dislocation–disclination field theory has been used to describe a hierarchy of cellular microstructures in ductile metallic crystals under severe plastic deformation [9, 53]. In this context, low-angle misorientations across cell walls at a relatively small length scale, termed incidental dislocation boundaries, are resolved by GNDs. Disclinations resolve potentially higher angle misorientations across larger cell blocks separated by geometrically necessary boundaries. The disclination concept thereby introduces another length scale for regularization into the theory and in conjunction with GNDs enables a natural description of microstructure evolution at different scales. The theory of [3, 53], thought to be the first thermodynamically complete constitutive framework for continuum dislocation–disclination mechanics at finite strain, includes a set of microscopic equilibrium equations, that in combination with kinematic relations between lattice spin and disclination density, enable prediction of the evolution of Q in an incremental boundary value problem. Numerical solutions of this formidable set of equations, which inevitably include entities from Riemannian and nonRiemannian geometry, remain elusive to date due to complexity. Furthermore, parameters are required that specify how free energy and slip resistance may depend on disclination density. Linear [3, 61] and nonlinear [64] elastic solutions for discrete defects offer insight into energies contributed by disclinations in certain arrangements, but unique assignment of parameters for hardening due to GNDs, SSDs, and disclinations remains problematic and highly material dependent. In contrast, more recent work has obtained numerical solutions for dislocation–disclination field theory applied to grain boundaries and triple junctions in metals [65], albeit with the constitutive theory restricted to the geometrically linear regime. A finite deformation description of the kinematics of GBs in relatively brittle minerals (i.e., low dislocation mobility) invoking a disclination density tensor has also been exercised [66]. A final concept is of importance when describing regions of very large defect density using continuum mechanical concepts. In such cases, the two term multiplicative decomposition of (2.9) may be
insufficient when FP is attributed to dislocation slip processes alone as implied in (2.10). An intermediate term, denoted here by FI , augments the classical crystal plasticity decomposition, accounting for residual lattice deformation due to defects within the local volume element of crystalline material at X to which F is ascribed: F ¼ r0 x ¼ FE FI FP :
ð2:26Þ
The particular form of FI depends on the class of defect, defect arrangement, and scale of resolution. A tensor equation for FI has been derived via homogenization (i.e., volume averaging) methods for polycrystals [29] and single crystals containing subgrain boundaries [59]. Elsewhere it has been calculated via consideration of the linear elastic fields of periodic arrays of edge dislocations [67, 68]. Solutions also exist for a volume element containing a single edge or screw dislocation [69, 70]. Representations in terms of nonlinear elasticity and anharmonic molecular statics have also been derived [3, 71], the former invoking third order elasticity of cubic crystals [72]. Analytical calculations have demonstrated the importance of inclusion of FI in the constitutive description when dislocation densities approach the theoretical maximum, an occurrence possible in regions of crystal near boundaries induced during severe plastic deformation or shock loading [70, 73]. In particular, since FP is isochoric when attributed solely to slip, any residual volume changes in the crystal are omitted if not captured by FI . Furthermore, lattice rotations induced by disclinations may be described by rotational part of FI [53], and volume changes associated with point defects may be included in its determinant [62, 74]. If included in a constitutive model framework, an evolution equation or one of the aforementioned subscale solutions for FI must be invoked. Though not essential, the residual lattice deformation may be included in the list of state variables fng, in which case additional model parameters may be needed to relate it to stored energy and strain hardening kinetics, for example. In summary, the classical sharp interface description invoked in crystal plasticity simulations, as presented in the ‘‘Sharp interface models’’ section, can be augmented to explicitly account for defect densities in the vicinity of GBs and subgrain boundaries. Potential benefits of the model classes covered here in the ‘‘Continuum dislocation and disclination
J Mater Sci
models’’ section include physical flexibility, where size effects and microstructure evolution at different resolutions can be incorporated somewhat naturally. Regularization associated with higher-order gradient terms in the free energy and balance or kinetic laws may facilitate mesh-size independence of numerical solutions, an issue which is of particular importance for modeling localization phenomena. A drawback is that additional parameters must be prescribed and often calibrated rather than determined from first principles, especially those relating defect densities to slip kinetics. Enhanced kinematics and stored energy, on the other hand, may be reasonably incorporated without calibrated parameters, via consideration of mathematical physics of defect densities (i.e., differential-geometric relations) and (non)linear elasticity solutions for individual non-interacting defects or those in idealized yet still sufficiently realistic arrangements. Regardless of the source of parameters, increased model sophistication results in an increase in computational expense and enables fewer, if any, available analytical solutions to boundary value problems for validation of computational results.
Diffuse interface models In the sharp interface classes of models in the ‘‘Sharp interface models’’ and ‘‘Continuum dislocation and disclination models’’ sections, lattice orientations demonstrate jump discontinuities across GBs, as do material properties such as elastic constants and thermal expansion coefficients that depend on the reference orientation of the crystal in Lagrangian constitutive models of anisotropic media. Even though regions of finite volume in the vicinity of boundary interfaces may contain distributions of defects represented by density tensors that may tend to smooth the mechanical response over small but finite distances from GBs, the corresponding models discussed in the ‘‘Continuum dislocation and disclination models’’ section still treat GB interfaces as discrete/sharp surfaces with regard to properties such as elastic coefficients. Diffuse interface models, perhaps most notably those termed phase field models, treat GBs as regions of finite volume over which referential properties vary continuously with distance from the interior of one crystal to the interior of its neighbor. Denote a scalar order parameter associated with a given grain
boundary shared by grains g1 and g2 by g 2 ½0; 1, such that gðX; tÞ ¼ 0 8X 2 Xðg1Þ ; gðX; tÞ 2 ð0; 1Þ
gðX; tÞ ¼ 1
8X 2 Xðg2Þ ;
8X 2 oXðg1;g2Þ ; ð2:27Þ
where now boundary zone oXðg1;g2Þ is of finite volume. Time is denoted by t. Then, a generic Lagrangian property A, which could be a tensor, vector, or scalar, is interpolated in GB regions from its constant values Aðg1Þ and Aðg2Þ initially assigned to regions deep within neighboring crystals as AðX; tÞ ¼ Aðg1Þ þ /½gðX; tÞsAt ¼ Aðg1Þ
ð2:28Þ
þ /½gðX; tÞðAðg2Þ Aðg1Þ Þ;
where / is an interpolation function minimally satisfying the end conditions /ð0Þ ¼ 0 and /ð1Þ ¼ 1. As discussed in the monograph [18] for example, an immense literature exists on diffuse interface models used to describe microstructure generation and evolution under various thermal and chemical processes: solidification, grain growth, recrystallization, mass transport, diffusion, etc. The scope of the present discussion in the ‘‘Diffuse interface models’’ section is hereafter limited to classes of models that treat GBs as diffuse but also address the mechanical response, specifically stress–strain behavior in crystals, which may be elastic or elastic–plastic. A relatively recent example of a coupled description of microstructure and mechanics via diffuse grain boundaries is reported in [75], which is focused on the electromechanical response of polycrystalline ferroelectrics. The phase field approach to lattice orientation assignment and GB representation follows that of [76], with essential equations of the method outlined in what follows next. A single order parameter suffices for description of one GB, e.g., that in a bicrystal. Diffuse interface models of GBs in polycrystals with many (n) grains require multiple order parameters, labeled gðiÞ , where i ¼ 1; 2; . . .; n. Following [75, 76], regions deep within each grain are characterized by values gðiÞ ðX; tÞ ¼ 1;
gðjÞ ðX; tÞ ¼ 0
8X 2 XðiÞ ;
ði 6¼ jÞ: ð2:29Þ
For any GB region, for at least one value of i,
J Mater Sci
jgðiÞ j 2 ð0; 1Þ:
ð2:30Þ
In the absence of mechanical (or thermal, electrical, etc.) loading, a free energy density per unit reference volume is prescribed as 1X fðgðiÞ ; r0 gðiÞ Þ ¼ f0 ðgðiÞ Þ þ ji jr0 gðiÞ j2 ; ð2:31Þ 2 i where ji are material constants. The smaller the value(s) of ji , the thinner the equilibrium width of a GB interfacial zone. The particular function f0 used in [75] is X a XX b f0 ¼ ðgðiÞ Þ2 þ ðgðiÞ Þ4 þ c ðgðiÞ Þ2 ðgðjÞ Þ2 ; 2 4 i i j6¼i ð2:32Þ with a, b, and c material properties. The function f0 contains 2n wells/minima at ðgð1Þ ; gð2Þ ; . . .; gðnÞ Þ ¼ ð1; 0; . . .; 0Þ, ð1; 0; . . .; 0Þ, ð0; 1; . . .; 0Þ; . . .. Function f0 is of a local minimum value within a grain, and r0 gðiÞ ¼ 0 8i within a grain, i.e., far from an interface. The free energy functional for the unloaded body is the integral over the entire polycrystalline volume X: Z F¼ fdX: ð2:33Þ X
Evolution of order parameters follows the Allen and Cahn [77] formalism, also known as the time-dependent Ginzburg–Landau (TDGL) equation, which drives the total free energy F of the system to a minimum. Letting l denote the mobility of the process (i.e., a parameter controlling the timescale of GB kinetics), the local rate equation for each order parameter in this approach is df of0 of ðiÞ g_ ¼ l ðiÞ ¼ l r0 dg or0 gðiÞ ogðiÞ 0 1 X ½gðjÞ 2 þ ji r20 gðiÞ A: ¼ l@agðiÞ b½gðiÞ 3 2cgðiÞ j6¼i
ð2:34Þ The following scalar function that has a value of unity in each grain and a magnitude less than unity within each GB region is introduced: X f½gðiÞ ðX; tÞ ¼ ½gðiÞ 2 2 ð0; 1: ð2:35Þ i
An interpolation function for properties used in [75] to define the fracture surface energy as GGB ¼ /GC of GBs is f itself:
/½fðgðiÞ ðX; tÞÞ ¼ fðgðiÞ ðX; tÞÞ 2 ð0; 1:
ð2:36Þ
This approach assigns the critical energy release rate property to a point X in a GB, denoted by GGB , as some positive fraction of the constant critical energy release rate GC assigned to all bulk crystals. Furthermore, in 2D simulations, an orientation angle H is defined at any point X via interpolation as [75] 1 X ðiÞ ðiÞ HðX; tÞ ¼ H ½g ðX; tÞ2 : ð2:37Þ fðX; tÞ i Orientation-dependent Lagrangian properties such as anisotropic elastic constants are then assigned to GB regions based on the local value of H, where HðiÞ is the uniform angular orientation of grain i. Notice that according to this model, lattice orientation in a given GB region potentially depends on orientations of (many other) grains not in contact with that boundary. The diffuse interface framework described in (2.29)–(2.37) is invoked in [75] to create a polycrystalline microstructure and assign orientations and fracture strengths to all points within the domain, i.e., GB regions and bulk crystalline regions. The TDGL equation is solved over a finite domain in time, beginning with randomly seeded admissible values of gðiÞ at points X 2 X used as initial conditions, from which grain growth take place. Time integration of (2.34) must be halted when a realistic microstructure (e.g., a realistic average grain size) is obtained; otherwise, a uniform single crystal would ultimately produce the minimum value of total system energy F. The polycrystal is then subjected to electromechanical loading, with the reference configuration (i.e., microstructure) held fixed in simulations. Another example of diffuse interface modeling of GBs in deformable polycrystals is reported in [78, 79]. Crystal elastic–plastic theory is used to represent the deformation behavior of aluminum grains in the aggregate. Accumulated dislocations and associated stored energy then supply driving forces for GB motion, e.g., grain growth during recrystallization. A staggered numerical scheme is implemented to solve the governing equations of crystal plasticity and phase field kinetics, where the solution of one set of physical laws influences that of the other set in successive iterations. Advantages of the diffuse interface models of GBs include the following. Microstructure evolution, e.g.,
J Mater Sci
GB motion, can be addressed in a more detailed and realistic way than in phenomenological models with sharp interfaces. In particular, kinetics of microstructure evolution are motivated by the fundamental principle that a system should seek a configuration for which its total free energy is a minimum. Regularization associated with gradient terms in the energy functional introduces length scale effects and facilitates mesh-size independence of numerical solutions. Disadvantages of the diffuse interface models include requisite prescription of non-unique interpolation functions for property values within interfacial zones such as / of (2.36), as well as parameters in the energy functional and kinetic law such as (ji ; a; b; c; l) in the framework of [75]. In computer simulations, mesh resolution must be fine enough to resolve field variables and their gradients within the often very narrow GB zones. Simultaneous solution of the coupled governing equations for the microstructure-mechanics problem may be challenging, especially if vastly different timescales for GB evolution and stress dynamics (e.g., wave propagation) arise.
Twinning Twins are microstructure features observed in many kinds of crystals. Twinning may be caused by mechanical forces, in which case it is termed deformation twinning or mechanical twinning. Twins may also be induced by other physical stimuli, a notable example being annealing twins produced via thermal processing of a material. The present section focuses mostly on deformation twins, particularly model descriptions of pseudo-slip and phase field type in ‘‘Continuum pseudo-slip models’’ and ‘‘Diffuse interface models’’ sections, respectively. Many aspects of sharp interface models discussed in the ‘‘Sharp interface models’’ section may apply to nearly any kind of twin, regardless of its origin. Prior to presentation and evaluation of the aforementioned classes of twinning models, a few fundamental concepts are reviewed. More complete treatments of twinning in the context of elastic–plastic continuum mechanics include [3, 10, 11]. Twinning is a general term that may be used to describe energy invariant transformations of a crystal structure with certain characteristics. A twin in a crystalline solid is usually defined as two regions of a
crystal separated by a coherent planar interface called a twin boundary. As will be described mathematically in the ‘‘Sharp interface models’’ section, limiting values of deformation gradients in each region, on either side of the twin boundary interface, differ by a simple shear. Unstressed twinned regions of the crystal far from boundaries or defects possess the same strain energy density as the unstressed parent (i.e., the same energy density as the original crystal prior to twinning), such that twinning shears are said to be energy invariant [80, 81]. In the context of ductile solids, deformation twinning is most often associated with thermodynamically irreversible shape deformation in correspondence with collective motion of partial dislocations and formation of stacking faults [11, 82]. Deformation twinning is preferred over slip in cases wherein resistances to dislocation glide are very large in certain directions, often in crystal systems of low, e.g., non-cubic, symmetry. In addition to their emergence in ductile metals, deformation twins may also appear in ceramics, minerals, and molecular crystals [83], though complex low-symmetry crystal structures do not ensure their occurrence [84]. Twinning is often preferable to slip at lower temperatures or at very high strain rates, though exceptions are not unusual, depending on material. Mechanical work done during deformation twinning is dissipative when resulting from defect motion associated with shearing. Any stored energy is associated only with defects left behind in the crystal, for example those comprising the twin boundary. From the standpoint of continuum thermodynamics, the driving force for twin propagation is the resolved shear stress on the habit plane in the direction of twinning shear, as will be exploited in the context of pseudo-slip models in the ‘‘Continuum pseudo-slip models’’ section.
Sharp interface models A sharp interface model of a twin boundary (TB) is similar to a sharp interface model of GB. In the latter, as discussed in the ‘‘Sharp interface models’’ section, a planar interface oXðg1;g2Þ separates two distinct crystals g1 and g2 with different referential lattice orientations, and Lagrangian material properties are discontinuous across the GB. In the former (TB), a planar interface oXðp;t1Þ separates the original crystal (parent) p and the twinned crystal t1, or more
J Mater Sci
generally, two different twins if p is replaced with t2. Lattice orientations demonstrate a jump discontinuity across the TB, as do associated anisotropic Lagrangian properties such as elastic moduli. However, the misorientation across a TB is restricted by the crystal structure and type of twin, while that across a GB is relatively unrestricted. Furthermore, a simple shearing process describes the transformation from the original lattice to the twinned lattice, whereas no such process generally exists for an arbitrary GB. Unlike general GBs that have a finite radius of curvature, fully formed TBs tend to be flat, though exceptions are common for growing or receding twins or those induced by concentrated forces [10] such as those encountered in (nano)indentation. The present focus is restricted to coherent TBs, for which continuity of displacement and traction hold analogously to (2.18) and (2.19): suðX; tÞt ¼ 0 , sxðX; tÞt ¼ 0; ð8X 2 XðpÞ ; Xðt1Þ and
stt ¼ srtn ¼ 0;
ð8Xðx; tÞ on oX
Þ:
W½Fðt1Þ ðXÞ ¼ W½Q0 FðpÞ ðXÞH;
ð3:5Þ
where Q0 is a proper orthogonal tensor T (Q1 0 ¼ Q0 ; det Q0 ¼ 1) and H is an energy invariant transformation of the crystal, not necessarily orthogonal, that depends on the material’s intrinsic structure/symmetry. In a single global Cartesian coordinate system, assuming that H does not induce volume changes since the converse assertion would permit the total strain energy of a sample of fixed mass to remain constant as the volume of the sample is increased without bound, and noting that det Fðt1Þ ¼ det FðpÞ ¼ det½FðpÞ þ a m0
ð3:2Þ
¼ det FðpÞ ½1 þ c0 s0 m0 :
parent p and twin t1, be separated by a surface oXðp;t1Þ across which displacements of the material are continuous as in (3.1), but across which gradients of displacement are not. This surface of composition [80], which need not be planar, corresponds to the habit plane in the traditional description of mechanical twinning [82]. Let Fðt1Þ ¼ r0 xðt1Þ
where c0 [ 0 is a scalar magnitude of the twinning deformation (eigen-shear) and s is a spatial unit vector. Let WðFÞ denote the strain energy density per unit reference volume of the crystal. Energy invariance of twinning demands that
ð3:1Þ
Kinematics of twinning can be described, in part, by invoking geometrically nonlinear elasticity theory. Let two regions of the crystal, which are labeled as
FðpÞ ¼ r0 xðpÞ ;
ð3:4Þ
det H ¼ 1, at the same limiting point X 2 oXðp;t1Þ ,
8X on oXðp;t1Þ Þ; ðp;t1Þ
sFt ¼ Fðt1Þ FðpÞ ¼ a m ¼ c0 s m0 ;
ð3:3Þ
denote constant limiting values of deformation gradient FðX; tÞ in each region in the vicinity of the TB, where X are reference coordinates of the original crystal prior to twinning. Since volumes and masses remain positive, det FðpÞ [ 0 and det Fðt1Þ [ 0. Let m0 be a unit normal vector to oXðp;t1Þ , pointing from parent side to the twinned side. The compatibility requirement that the interface be coherent (i.e., continuous coordinates x along the surface of composition or TB) necessitates that Hadamard’s jump conditions apply [1, 80, 85]:
ð3:6Þ
It follows that the pullback of s, i.e., s0 ¼ ½FðpÞ 1 s, must be orthogonal to unit normal m0 : s0 m0 ¼ 0:
ð3:7Þ
When the parent is taken as a perfect reference lattice, then FðpÞ ¼ 1 ) Fðt1Þ ¼ 1 þ c0 s0 m0 ;
ð3:8Þ
demonstrating that Fðt1Þ is indeed a simple shear. In that case, the equivalent product Q0 H is also a simple shear, possibly of large magnitude, that shifts the perfect crystal to another minimum energy configuration, with the strain energy density of this configuration equivalent to that of the parent. In this context, the strain energy function W can be interpreted as a multi-well potential, with global minima corresponding to conditions Wð1Þ ¼ WðQ0 HÞ ¼ 0, a characteristic feature that will be exploited later in the ‘‘Diffuse interface models’’ section in the context of diffuse interface theories of twinning. Notice that the above description does not account for any (surface) energy associated with defects along the boundary of the twin, which can be reflected in continuum theories via augmentation of the free energy function with internal state variables, as will be demonstrated in the ‘‘Continuum pseudo-slip models’’ section.
J Mater Sci
The preceding treatment addresses kinematics and strain energy density for sharp interface model representations of twin boundaries. Elements of such a treatment can be invoked to analyze and predict occurrence of various preferred microstructures. For example, laminated twin arrangements (i.e., herringbone patterns) and other characteristic features of martensite [1, 85] can be predicted from consideration of compatibility constraints and free energy minimization for certain crystal structures. No material parameters generally need to be calibrated in such analyses, which rely on fundamental properties such as symmetry operations—and transformation strains if solid–solid phase changes are involved—associated with crystal structure. The above treatment does not enable explicit prediction of time-dependent motion of twin boundaries, e.g., twin growth and dynamic interactions with other twins or other crystals in a polycrystal. For such predictions, a kinetic law for twin boundary dynamics must supplement the kinematic description, and a (numerical) scheme must be invoked to track the position of interface(s) as deformation proceeds in time. Derivation of the corresponding equations is beyond the present scope, but one such example of this class of sharp interface model for twinning dynamics is the 2-D theory and level-set numerical method of [86]. Therein, a stored energy function is non-convex with multiple wells. Evolution of twin interfaces is governed by a kinetic relation for the twin boundary velocity as a function of the local driving traction and boundary orientation. A regularized version of the theory is constructed via the level-set method which removes the requirement of treatment of explicit jump conditions. Numerical finite difference results in [86] compare favorably with observed phenomena in martensite: cusp formation, needle growth, spontaneous tip splitting, and microstructure refinement. Unlike purely (nonlinear) elastic treatments, however, sharp interface models of TB motion require kinetic laws and material parameter(s) that relate driving forces to interface velocities, for example.
Continuum pseudo-slip models In what is termed here as a pseudo-slip class of model, a local volume element of a crystal consists of fractions of the parent and one or more deformation twins. Volume or mass fractions of twins evolve
according to a kinetic law, with the driving force for twinning typically a resolved shear stress acting on the habit plane of the twin system, in the direction of twinning shear for that system. Twin boundary interfaces are not resolved explicitly within each volume element. However, the boundary between a fully twinned domain and the parent or a domain containing twins of other twin systems is captured in a homogenized or smoothed sense, whereby neighboring coordinates X may support different volume fractions of each twin variant. The pseudo-slip approach was apparently first introduced in [87, 88] where it was used for crystallographic texture predictions. Finite element implementations of a purely mechanical theory accounting for elasticity, slip, and twinning were perhaps first reported in [89, 90], where pseudo-slip laws were invoked for ductile metals. The first complete thermomechanical frameworks accounting for such deformation phenomena—both exercised to describe shock or high-pressure phenomena for which higherorder thermoelasticity is essential—are described in [73, 83]. The first to also include GNDs in a gradient theory, merging the nonlinear thermoelasticity, crystal plasticity, and pseudo-slip twinning descriptions with concepts described in ‘‘Continuum dislocation and disclination models’’ section, was presented in [91]. Results of these works address ceramics and minerals [73, 91] or molecular energetic crystals [83]. The forthcoming presentation summarizes the theory developed in [3, 73, 91]. The deformation gradient is decomposed into a product of three terms: F ¼ r0 x ¼ FE Fg FP ;
ð3:9Þ
where thermoelastic deformation FE and deformation from plastic slip FP have the same meanings as in (2.9) of the ‘‘Continuum dislocation and disclination models’’ section. The contribution of twinning shear to the total deformation gradient for a volume element at point X and time t is denoted here by Fg ðX; tÞ. A term akin to det FI of (2.26) is also used in the full kinematic framework of [3, 73, 91] to account for possible volume changes associated with lattice defects (including dislocation cores, stacking faults, and TBs) but is omitted here in the interest of brevity. Let gb ðX; tÞ denote the volume fraction of twin variant b of the material at point X and time t, where b ¼ 1; 2; . . .; q, with q the number of twin systems. Let cb0 denote the stress free twinning shear associated
J Mater Sci
with variant b, which has shearing direction sb0 and mb0 .
Then, the rate of deformation habit plane normal and spin from twinning is computed via X Lg ¼ F_ g Fg1 ¼ g_ b cb0 sb0 mb0 : ð3:10Þ b
The plastic velocity gradient of (2.10) is replaced with an augmented equation that accounts for the change in lattice orientation of twinned domains: LP ¼ Fg F_ P FP1 Fg1 ¼ ð1 gT Þ
X
c_a sa ma þ
a
X b
gb
X
! c_ab sab mab :
a
ð3:11Þ P
gb 2 ½0; 1 is the total twinned volume Here, gT ¼ fraction, and the second (double) sum is over twinned domains with slip rates c_ab and rotated or reflected director vectors sab and mab . The rotation or reflection matrices depend on the crystal structure and twin type, e.g., a type I or type II twin [3, 11]. Equations for thermoelasticity such as (2.11)–(2.13) still hold, but with the anisotropic thermoelastic moduli updated according to the weighted rotation/ reflection matrices to account for the twinning transformation. The free energy function is of the same generic form as in (2.14), but the twinning shears must be included in the list of internal state variables fng to enable description of the effect of twinning on the anisotropic thermoelastic coefficients. In the theory of [91], for example, fng ¼ fgb ; qG ; qS g:
ð3:12Þ
The evolution of the twin variants is dictated by a pseudo-slip law: g_ b ¼ g_ b ðsb ; T; fngÞ;
ð3:13Þ
where the twinning direction (sign of resolved shear stress sb acting on the variant’s habit plane) must be respected to account for increased deformation resistance in the anti-twinning sense. The local energy balance (i.e., temperature rate equation) accounts for dissipation from twinning shear in addition to that from plastic work, extending the elastic–plastic representation in (2.16). The pseudo-slip-based class of models to which the above theory and those developed in [83, 89, 90] belong enables reasonably accurate predictions of the onset and evolution of bulk twinning behavior and associated
crystallographic texture changes. Such models can be implemented in existing crystal plasticity simulation frameworks with modest additional effort. However, specific kinetic equations and parameters must be assigned to (3.13), and effects of twinning on slip must be included in the flow rule for the slip rates via augmentation of (2.15). Calibration and validation of such features are often problematic, with unique property selection difficult, if not impossible, due to the immense number of possible slip–slip, slip–twin, and twin–twin system interactions, each of which may most generally demonstrate different physical behaviors [11, 89]. Even more complexity is introduced if de-twinning is incorporated. Phenomenological expressions may be assigned to describe an evolving thickness to each local twin variant as in [91], but shapes of each twin variant are not predicted within a volume element at X by this class of models. If the gradient aspect of the theory (i.e., qG ) is omitted as in [73, 83, 89, 90], the model contains no intrinsic length for regularization; thus, those models’ predictions do not depend on the absolute size of the domain. Twin boundary migration has recently been explicitly incorporated in a computational crystal plasticity framework applied to nano-twinned metals [92].
Diffuse interface models The diffuse interface representation of a twin boundary (TB) has many similarities to the diffuse interface modeling scheme of GBs outlined in the ‘‘Diffuse interface models’’ section. One or more order parameter(s) are introduced that delineate the parent crystal from one or more twin variant(s). Deep within the parent and deep within each twin, order parameters are homogeneous, typically with numerical values of zero or unity depending on details of the model formulation. Twin boundaries are represented by finite volumes within which spatial gradients of order parameter(s) do not vanish. Physical properties that depend on lattice orientation (e.g., anisotropic Lagrangian elastic constants) are interpolated between parent and twin(s) across the boundary regions. Unlike the GB models in which deformation mechanics are not often addressed, in diffuse interface models of deformation twinning, accounting for the kinematics is paramount. Specifically important are the kinematics of the formation of the interface and the transformation (i.e., stress free shearing) of the twin variant(s).
J Mater Sci
The first nonlinear phase field theory for twinning in crystals, incorporating both nonlinear anisotropic elasticity and geometric nonlinearity, appears to be that of [93]. Around the same time [94, 95] appeared, albeit limited to small deformations, e.g., linear elasticity. In general applications, incorporation of nonlinear elasticity is deemed crucial, since different predictions arise in analytical [96] and numerical results [97, 98], and since shears, rotations, and/or reflections inherent to the twinning process all tend to be large, exceeding the usual limits of continuum linear elastic constitutive models. Supporting the asserted general necessity of nonlinear theory, Fig. 4 demonstrates different results for twinning in calcite single crystals modeled via the nonlinear theory of [93, 97] and its linearization. The lamellar features observed in many instances (e.g., in martensite) are present in the nonlinear result but absent in the linear result, speculatively due to some greater departure from convexity of the total potential energy in the former. Another nonlinear theory appearing soon after [93, 97] is applied to martensite in [99]. In what follows next, key features of the variational, finite deformation, phase field theory for deformation twinning developed and refined in [93, 97, 100] are reviewed. Attention here is limited to a single crystal with a single twin variant; extension to multiple twin systems is conceptually straightforward and is described in an appendix of [93]. This particular theory has not been used in conjunction with a model component for plastic slip (i.e., glide of dislocations distinct from twinning partials), meaning the response is limited to combined elastic
deformation and twinning deformation. However, an example of a coupled phase field-crystal plasticity model for metals that undergo simultaneous dislocation slip and twinning has been reported [101]. Thermal effects are omitted. Let the reference volume of a crystal X (which may be contained within a polycrystal) be divided into parent XðpÞ , twin Xðt1Þ , and boundary oXðp;t1Þ regions. The order parameter associated with twinning is denoted by g and obeys gðX; tÞ ¼ 0 8X 2 XðpÞ ; gðX; tÞ 2 ð0; 1Þ
gðX; tÞ ¼ 1
8X 2 Xðt1Þ ;
8X 2 oXðp;t1Þ : ð3:14Þ
Deformation and displacement are defined as in (2.8), and the continuity requirements in (3.1) and (3.2) still apply. However, unlike the sharp interface treatment of the ‘‘Sharp interface models’’ section, the deformation gradient now is presumed continuous everywhere, including points in oXðp;t1Þ . The total deformation gradient obeys decomposition F ¼ r0 x ¼ FE Fg ;
the ð3:15Þ
where FE is the elastic part and Fg accounts for shearing due to mechanical twinning. Let s0 and m0 denote constant orthogonal unit vectors in the direction of twinning shear and normal to the habit plane, respectively. Let c0 [ 0 denote the magnitude of stress free shear for a fully transformed domain. Then, Fg ðgÞ ¼ 1 þ ½/ðgÞc0 s0 m0 :
ð3:16Þ
The interpolation function /ðgÞ obeys /ðgÞ 2 ½0; 1; /ð0Þ ¼ 0; d/ d/ ð0Þ ¼ ð1Þ ¼ 0: dg dg
/ð1Þ ¼ 1; ð3:17Þ
Let w denote the free energy density per unit reference volume of the following form: wðF; g; r0 gÞ ¼ W½FE ðF; gÞ þ fðg; r0 gÞ;
ð3:18Þ
with W the elastic strain energy density. The function f, which is nonzero only in TB regions, consists of the sum of a double-well potential f0 and a gradient contribution: Figure 4 Phase field order parameter g representing twinning in nano-indentation of calcite single crystal with 120 wedge: nonlinear theory (left) and linear theory (right) [97].
J Mater Sci
fðg; r0 gÞ ¼ f0 ðgÞ þ j : r0 g r0 g 2
2
¼ Ag ð1 gÞ þ j : r0 g r0 g:
ð3:19Þ
Here, constant A quantifies the depth of the energy wells, and second-order tensor j penalizes sharp interfaces. When j ¼ j1, with j a scalar, the TB energy is isotropic. In this case, the equilibrium thickness l and equilibrium surface energy per unit reference area C are related to parameters in (3.19) by [93] A ¼ 12C=l;
j ¼ 3Cl=4:
ð3:20Þ
Let the total free energy functional be denoted by W, and let t0 denote the mechanical traction vector per unit reference area and h a conjugate force to the order parameter acting on global external boundary oX that has unit outward normal vector n. The following variational principle is applied: Z I I dW ¼ d wdX ¼ t0 dudoX þ hdgdoX: oX
X
oX
ð3:21Þ Application of standard mathematical techniques for analysis of continuous media then results in local equilibrium equations in X and natural boundary conditions on oX: r0 P ¼ r0
t0 ¼ P n;
oW ¼ 0; oF
df0 oW ¼ 2r0 jr0 g; þ og dg ð3:22Þ
h ¼ 2j : r0 g n:
ð3:23Þ
The first Piola–Kirchhoff stress is P ¼ ðdet FÞrFT . The constitutive model is complete upon specification of the strain energy function W and the interpolation function /. Regarding the former, any hyperelastic potential suitable for crystalline media can be used, but the elasticity tensor(s) C (if anisotropic) must be interpolated in the TB regions, e.g., C½gðX; tÞ; X ¼ CðpÞ ðXÞ þ /½gðX; tÞfCðt1Þ ðXÞ CðpÞ ðXÞg: ð3:24Þ Two interpolation functions have used with success to date, a cubic polynomial and a Fermi–Dirac exponential [98]: / ¼ ð3 2gÞg2 ;
/¼
1 : 1 þ exp½30ðg 1=2Þ ð3:25Þ
The latter function in (3.25) yields a rather steep change in / near g ¼ 1=2. Although this model is variational and quasi-static, evolution of twin morphology during deformation paths is predicted by sequential energy minimization (i.e., minimization of W subject to any essential boundary constraints) as loads are incrementally applied. The predictive capability of the theory described above has been validated or favorably compared with results from experiments, analysis, or atomic simulations for a number of crystalline materials—calcite, sapphire, and magnesium—for problems involving twin nucleation and/or growth from a seedling [93], indentation [97], and a notch or crack tip [98]. A generic twinning criteria based on analytical solutions to the phase field theory for localized versus diffuse transformation has also been validated [100]. The primary purpose of the diffuse interface approach to modeling twins and TBs is prediction of detailed, fine-scale twin morphology. Balance laws (or kinetic laws if a TDGL or Allen-Cahn [77] equation is invoked as in [94, 95]) are derived from fundamental principles of energy minimization rather than user-prescribed phenomenology as is typical in pseudo-slip models of the ‘‘Continuum pseudo-slip models’’ section. Furthermore, few, if any, material parameters must be calibrated. This advantage partially disappears if continuum slip laws requiring parameterization are added to the theory, as in [101], for example. The length scale of regularization l is controlled by the modeler and dictates the equilibrium width of TB zones. The regularization process renders numerical simulations mesh-size independent so long as the discretization is fine enough to resolve order parameter gradients in such zones. This is often a drawback in simulations involving twin(s) of size much larger than their boundary zones, since very fine meshes must be used to resolve boundaries that encompass a relatively small overall fraction of the entire problem domain. The motion of twin or crystallographic phase boundaries—be it via dynamic extension, thickening, or migration—may crucially influence the mechanical response of certain crystals, particularly nano-twinned metals [102], shape memory alloys, and martensite [85, 103]. A novel theory and complementary computational method were recently developed that includes distinct prescriptions of
J Mater Sci
nucleation through the source term of the phase field conservation law and kinetics through a distinct interfacial velocity field [104, 105]. This approach alleviates potential obscurity of these physical phenomena suffered by classical phase field implementations.
Failure and localization Mesoscopic continuum models of damage processes in polycrystals are now addressed. These processes tend to be irreversible and dissipative. Examples include fracture, void formation and growth, pore collapse, and shear localization. In crystalline materials, fractures may occur along GBs, i.e., intergranular modes, and/or along cleavage planes within a grain, i.e., transgranular modes. Regarding intergranular fracture, some kinds of crystals demonstrate preferred planes, typically of low surface energy, while others do not, in which case fractures tend to be conchoidal [106]. Void mechanisms arise in ductile as opposed to brittle solids, wherein plastic deformation of surrounding matrix material enables volumetric growth of damage rather than planar fractures. In materials with initial porosity, including many minerals and their composites (e.g., concrete [107]) pores may be irreversibly compressed out of the material due to mechanical pressure, leading to an increase in compressive bulk modulus. As discussed already in the ‘‘Sharp interface models’’ section, shear bands may arise due to lattice orientation effects [39], but are perhaps most common in high rate deformation processes wherein localization is promoted by (nearly) adiabatic heating and thermal softening [8]. Tensile wave interactions also can induce damage mechanisms, specifically spall failure [108], which may entail dynamic brittle fracture or ductile failure via void coalescence, depending on ductility of the polycrystal.
Cohesive fracture models Cohesive fracture models are a kind of sharp interface representation of failure at the mesoscale. The theoretical concept of such models, most often attributed to [109, 110], is that along surfaces near a crack tip, the degraded material supports a nonzero traction vector over some finite distance, called the cohesive zone length. The traction tends to decrease
in magnitude as the crack opening displacement increases. Traction generally includes both normal and shear components, in association with mode I and mode II/III crack opening dispacements [7]. The earliest numerical implementation of a cohesive failure model in a finite element context appears to be documented in [111]. Subsequently, cohesive failure elements have been invoked in quasi-static crystal plasticity simulations [112, 113], thermoelastic simulations of ceramics [44], and elastodynamic simulations of crack branching [114] and spall [45]. The first dynamic simulations coupling finite crystal thermoelasto-plasticity with cohesive failure along GBs seem to be those reported in [115, 116], with a followup study of spall in [117]. Representative simulation results for cohesive failure modeling of dynamic fracture in ceramic and metallic polycrystalline microstructures are shown in Figs. 5 and 6, respectively. In all such simulations, the bulk material response within continuum finite elements is modeled via standard thermoelasticity or crystal elastoplasticity, e.g, models discussed in the ‘‘Grain boundaries’’ section of this work. The cohesive zone model is invoked for failure/separation behavior along element boundaries. A few key equations for a basic cohesive zone model of fracture are reviewed next. An immense literature on cohesive fracture modeling has emerged over the previous two decades, with more sophisticated theoretical and numerical formulations accounting for various subscale physical mechanisms and thermodynamic aspects now available. The presentation below is the minimum deemed necessary to illustrate the fundamental mechanical concepts. Let the crack opening displacement vector across two crack faces initially coincident at point X be defined as the displacement jump dðX; tÞ ¼ suðX; tÞt;
ð4:1Þ
where now obviously the continuity constraint in (2.18) is violated across the crack surfaces. Denote by nðXÞ the unit normal to the surface of impending fracture, referred to the reference coordinate system. Let t0 denote the traction vector per unit reference area: t0 ¼ P n;
ð4:2Þ
with P the first Piola–Kirchhoff stress tensor. In the cohesive zone, a general form of traction–separation law is prescribed:
J Mater Sci
Figure 5 Shear stress r in AlON ceramic under simulated shear and compression with nonlinear elastic grains and cohesive finite elements at GBs: external view (left) and internal view (right) [27].
Figure 6 Particle velocity up in shock wave prior to spall (left, arrows are velocity vectors; colors indicate velocity magnitude) and effective stress at spall zones (right) in dynamic crystal plasticity-cohesive finite element simulations [117].
t0 ¼ t0 ðd; fvgÞ:
ð4:3Þ
Potential history effects are captured by state variable(s) fvg. Typically, a magnitude of displacement dC is assigned as a material property or parameter, beyond which traction vanishes and the crack surfaces become free surfaces. For quasi-static implementations, a stiffness matrix is usually needed, in which case the traction function must have a continuous derivative with respect to opening displacement. For simulations invoking explicit numerical integration, on the other hand, this differentiability restriction does not hold. The work done during separation can be related to a surface energy of fracture C, equal to half the critical strain energy release rate in the context of linear elastic fracture mechanics:
1 C¼ 2
Z
dC
t0 dd;
ð4:4Þ
0
where the path of integration ends when the critical separation magnitude is attained. Perhaps the simplest realistic model is the triangular cohesive law jt0 j ¼ rC ð1 jdj=dC Þ;
ð4:5Þ
which can be invoked separately for magnitudes j j of normal and shear components. Here, rC is the strength required to initiate fracture, i.e., the resolved scalar stress component at which the cohesive zone starts to open. Only two of the three parameters (C; rC ; dC ) need be prescribed since (4.4) enables one of these to be eliminated algebraically. More sophisticated models accounting for mode mixity are typical [114, 118], but these often need additional calibration or parameters. The length of the cohesive zone in the context of isotropic linear elastic fracture mechanics is [7, 118]
J Mater Sci
lC
2EC pr2C ð1
mÞ
2
;
ð4:6Þ
where E and m are the elastic modulus and Poisson’s ratio. Cohesive failure models often enable realistic predictions of brittle or ductile failure in microstructures. Distinct fracture entities, i.e., crack sizes and shapes, are fully resolved, and interactions among entities are naturally addressed. Advantages of this class of model include relatively few parameters [minimally two, e.g., dC and rC in the context of (4.5)] and distinct behavior of interfaces and bulk material, meaning that traditional solid continuum elements can be used for representing the latter. However, real polycrystalline solids demonstrate a distribution of strengths and surface energies among potential failure sites (e.g., various kinds of GBs, different families of cleavage planes, and initial defect structures that affect toughness), a characteristic often ignored in deterministic computer simulations. Numerical implementation of the basic model is rather straightforward, though additional nodal degrees of freedom increase the computational expense where duplicate nodes are inserted along failure planes. A potentially severe drawback is that fracture paths are constrained to follow element boundaries in most computer implementations. Therefore, meshes must be constructed such that realistic crack paths are possible, often requiring some a priori knowledge of failure morphologies. This drawback is alleviated when fractures are restricted to preexisting interfaces such as GBs, which makes the approach ideal for representing intergranular failure or fracture between phases of heterogeneous crystals such as reported in [115–117]. Cleavage fracture on specific planes is more difficult to address via cohesive element modeling, but examples of success have been reported [119]. Finally, cohesive finite element sizes must be small enough to resolve behavior over the cohesive length, e.g., must be smaller than lC of (4.6). Extremely fine meshes are often necessary for representing fracture zones of materials with high strength and low surface energy, e.g., many strong yet brittle ceramics, as modeled, for example, in [120].
Continuum damage mechanics In continuum damage mechanics models, individual failure entities such as discrete cracks and voids are
not resolved explicitly. Instead, one or more state variable(s) is introduced as a function of time t and material coordinates X that accounts for degradation of the material at a local material point and a given time instant. Such state variables, which may be scalar-, vector-, or tensor-valued, quantify the effects of multiple damage entities contained within a local volume element of material at that point. Research books on the subject include [121, 122]. The fundamental concept, in the context of a scalar damage mechanics theory, is often attributed to [123]. The forthcoming presentation considers a (poly)crystalline solid whose damage is represented by a scalar state variable DðX; tÞ 2 ½0; 1. Generalization to a vector, tensor, or multiple scalars is in principle straightforward but notationally more cumbersome. The solid is assumed to be hyperelastic and may undergo plastic slip, following the finite deformation continuum formalism discussed in the ‘‘Sharp interface models’’ and ‘‘Continuum dislocation and disclination models’’ sections. The deformation gradient is decomposed multiplicatively as F ¼ r0 x ¼ FE FD FP ;
ð4:7Þ
where FE includes thermoelastic deformation of the crystal lattice as well as mechanically reversible changes in damage (e.g., elastic crack closure on load release), FP accounts for plastic slip from dislocations as described in the ‘‘Sharp interface models’’ section, and FD accounts for mechanically irreversible damage mechanisms such as cracks and voids that remain in the solid upon elastic unloading. A three term decomposition of this general form was proposed in [124]. Other forms of the deformation gradient that reflect residual damage modes include additive [115, 125, 126] and hybrid additive–multiplicative [113, 127] decompositions, often derived or motivated from homogenization of discrete displacement jumps associated with cracks within a volume element of material. Transformations between additive and multiplicative descriptions have also been derived [125, 128]. The volume fraction of damage N (i.e., porosity) is related to the determinant of FD as [3, 129] J D ¼ det FD ¼ ð1 NÞ1=3 :
ð4:8Þ
Besides its use for solids with voids [124, 129] or pores [21, 107], a multiplicative damage term has been introduced for cleavage cracking in crystals
J Mater Sci
driven by resolved normal and shearing tractions in pseudo-slip-type models [32, 128, 130]. The list of internal state variables of (2.22) is extended to include damage in addition to dislocations: fng ¼ fqG ; qS ; Dg:
ð4:9Þ
Damage variable D varies from zero to unity as the material at the corresponding point loses integrity. The free energy density w therefore depends on damage in addition to elastic strain E, temperature T, and dislocation densities. The elastic second Piola– Kirchhoff stress and conjugate force to damage are S¼
ow ; oE
f¼
ow : oD
ð4:10Þ
An equation similar in form to (2.13) applies for the stress, but now elastic moduli depend on damage. The simplest degradation model of the moduli is linear in D: C½DðX; tÞ; X ¼ ½1 DðX; tÞC0 ðXÞ;
ð4:11Þ
with C0 ðXÞ ¼ Cð0; XÞ the tensor of elastic moduli for the undamaged crystal at the corresponding material point. For an isotropic solid, (4.11) leads to a constant Poisson’s ratio. More sophisticated degradation laws are required to more realistically capture physics of arbitrary loading paths, such as damage induced anisotropy and differences in tensile versus compressive degradation, where generally the former is more severe. The work done by the appropriate stress tensor acting on the rate of FD contributes to local dissipa_ Kinetic equations must tion, as does the product fD. be supplied for time rates of FD and D, e.g., of general forms F_ D ¼ F_ D ðE; fng; T; FD Þ;
_ D_ ¼ DðE; fng; T; FD Þ: ð4:12Þ
Dependence on elastic strain and internal state variables is often conveniently replaced with physically more transparent dependence on stress and other thermodynamic driving forces. Kinetic equations for plasticity and dislocation density are likewise affected by nonzero damage [124, 130]; e.g., stress concentrations may increase the tendency for plastic flow. Capabilities and caveats of the continuum damage classes of material models can be summarized as follows. The present class of models to which the
above theory, and typical of those discussed in [121, 122], enables reasonably accurate predictions of the onset and evolution of bulk damage behavior and associated changes in elastic stiffness. These models can often be implemented in existing continuum mechanics simulation frameworks, including those accounting for crystal plasticity, with modest overhead. Unlike cohesive finite element approaches, no special interfacial elements or node duplications are necessary. Unfortunately, specific kinetic equations and parameters must be assigned to (4.12), and coupled effects of D must be included in the flow rule for the slip rates via augmentation of (2.15). Calibration and validation of these aspects of the model are often difficult, and unique specification of all parameters cannot usually be ensured from available test data for a particular material. Additional complexity is injected when more realistic formulations involving multiple damage variables or vector- or tensor-valued damage variables are employed for anisotropic media. Since damage is homogenized at each point X, this class of models is unable to resolve solution fields of discrete cracks or voids and their interactions. The damage components of this class of model contain no intrinsic length for regularization, and problems that involve the usual material softening with increasing D often suffer from mesh-size-dependent numerical solutions. Such problems can be alleviated when gradient terms are considered, e.g., qG if GNDs are addressed as in (4.9) or [130]. Gradients of the damage variable can also be introduced into the balance equations and/or kinetics, which renders the continuum damage theory nonlocal and hence non-classical. This type of regularization is a characteristic advantage of diffuse interface or phase field approaches discussed next in ‘‘Diffuse interface models’’ section.
Diffuse interface models Diffuse interface classes of models for failure mechanisms, which include phase field representations of fracture, do not resolve discrete jumps in the displacement field in contrast to sharp interface-type (e.g., cohesive zone) models. Instead, fracture surfaces (or other failure surfaces) are described by one or more order parameter(s). In the usual case, an order parameter is assigned a uniform value (e.g., zero or unity) in regions of material where no damage exists or where complete failure has occurred.
J Mater Sci
The change from initial or undamaged state to fully failed state can be thought of as a transition from a perfect crystal to a liquid state or a vacuum, depending on the strength/stiffness properties attributed to the fully failed state. Boundaries between fully failed and perfectly intact material are represented by order parameter gradients, and in such boundary regions the order parameter takes on a value between its two extremes. Lagrangian elastic properties depend on the order parameter, with a decrease in local tangent stiffness correlating to an increase in local damage. Two of the earliest diffuse interface models of fracture are reported in [131, 132]. These early models invoked linearity in terms of deformation kinematics and elastic response. A nonlinear diffuse interface model for isotropic solids is presented in [133] in concordance with 2-D simulations. The first nonlinear phase field model for fracture in anisotropic crystals is described in [134]. In that work, 3-D simulations for mode I and mode II loading are validated versus analytical solutions, and a problem associated with crack bridging is investigated. Specifically, the tendency for a cleavage fracture to propagate around a spherical inclusion or cut through the inclusion is quantified in terms of stiffness and strength properties of the matrix and the sphere. A particular result
in which the crack is deflected around the inclusion is shown in Fig. 7. Mesoscale phase field simulations of simultaneous intergranular and transgranular fractures in anisotropic polycrystals, with and without secondary GB phases, are reported in [30], with characteristic results shown in Fig. 8. The forthcoming discussion presents key aspects of the nonlinear anisotropic phase field model for fracture of [134]. As noted above, displacement and traction are continuous fields, so constraints (2.18) and (2.19) hold. Plastic deformation, twinning, and thermal effects are omitted in what follows, noting that a coupled nonlinear phase field theory for simultaneous fracture and twinning mechanisms has been developed and analyzed in [100], with numerical simulation results for polycrystals given in [135]. Let the reference volume of a crystal X be divided into perfect XðpÞ , fully failed XðfÞ , and boundary oXðp;fÞ regions. The theory considers a single order parameter gðX; tÞ, where now gðX; tÞ ¼ 0 8X 2 XðpÞ ; gðX; tÞ 2 ð0; 1Þ
gðX; tÞ ¼ 1
8X 2 XðfÞ ;
8X 2 oXðp;fÞ : ð4:13Þ
Deformation and displacement are defined as in (2.8), with the deformation gradient and its determinant associated with volume change given by the usual equations F ¼ r0 x;
J ¼ det F [ 0:
ð4:14Þ
The free energy density per unit reference volume is of the following form: wðF; g; r0 gÞ ¼ WðF; gÞ þ fðg; r0 gÞ;
ð4:15Þ
with W the elastic strain energy density. The function f, which vanishes identically in the perfect/pristine regions of the crystal(s), consists of the sum of a quadratic potential f0 and a quadratic gradient contribution: fðg; r0 gÞ ¼ f0 ðgÞ þ j : r0 g r0 g C ¼ g2 þ fCl½1 þ bð1 m mÞg : r0 g r0 g: l ð4:16Þ
Figure 7 Phase field prediction of crack deflection around a strong second-phase inclusion for far-field mode I loading with failed material (g [ 0:7) removed to visualize crack propagation [134].
Here, C is the surface energy of fracture, l is the equilibrium width of the diffuse fracture zone, m is a unit normal vector to a preferred cleavage plane in the reference configuration, and b is a penalty factor that should be prescribed as a large number to restrict
J Mater Sci
Figure 8 Phase field predictions of inter- and transgranular fracture in anisotropic Zn polycrystal of edge length L under axial tensile strain e with failed material (g [ 0:7) removed to visualize
crack propagation [30]: a L ¼ 100 lm, e ¼ 0:15% b L ¼ 100 lm, e ¼ 0:20% c L ¼ 100 nm, e ¼ 4:4% d L ¼ 100 nm, e ¼ 8:0%.
fractures to take place along m. Isotropic fracture corresponds to b ¼ 0, in which case m is not needed. The strain energy density W may be any nonlinear elastic potential suitably modified to account for degradation of stiffness with increasing values of g. For isotropic neo-Hookean elasticity considered in some simulations [134], the shear modulus l degrades from its initial value l0 as
mechanical traction vector per unit reference area, and h is a conjugate force to the order parameter acting on boundary oX that has unit outward normal vector n. The first Piola–Kirchhoff stress tensor is P. A variational principle is applied: Z I I dW ¼ d wdX ¼ t0 dudoX þ hdgdoX;
lðgÞ ¼ l0 ½f þ ð1 fÞð1 gÞ2 ;
ð0 f 1Þ: ð4:17Þ
Parameter f, if nonzero, provides for finite stiffness in fully fractured zones. The bulk modulus degrades in tension but not in compression via a criterion based on local values of volume ratio J. Similar treatments for degradation of anisotropic elastic constants for crystals are described in [30]. Derivation of the Euler–Lagrange equations parallels that of the ‘‘Diffuse interface models’’ section. The total free energy functional is W, t0 denotes the
X
oX
oX
ð4:18Þ which results in local equilibrium equations and natural boundary conditions: r0 P ¼ r0
t0 ¼ P n;
oW ¼ 0; oF
2
Cg oW þ ¼ 2r0 jr0 g; l og ð4:19Þ
h ¼ 2j : r0 g n:
ð4:20Þ
Diffuse interface models of fracture or material failure exhibit some important positive characteristics: structural transformations occur naturally via energy minimization as in the approach outlined
J Mater Sci
above, or Ginzburg–Landau-type kinetics [131], without recourse to phenomenology associated with often ad hoc, user-prescribed rate equations for damage evolution. Correspondingly, relatively few material parameters or properties are usually needed. For example, in the theory of [30, 134] described above, the only required parameters are the elastic constants, fracture surface energy C, and the length constant l. Regarding the latter, numerical solutions are regularized and thus are rendered mesh-size independent by l, which is associated with interfacial width. Regularized variational models invoking the form of modulus degradation in (4.17) may also demonstrate so-called gamma convergence toward the classical Griffith theory of fracture mechanics as the regularization zone shrinks to a singular surface. Detailed damage morphologies can be predicted by numerical simulations invoking phase field fracture theory, including crack nucleation and branching and interactions among multiple cracks [136]. Standard continuum finite element technologies and meshing strategies can be used since fields are continuous, but additional degree(s) of freedom [i.e., order parameter(s)] must be updated at each node as the calculation proceeds incrementally. A disadvantage of the phase field approach to modeling fracture is the requirement of resolution of damage surfaces/ boundaries. A very fine mesh is needed to resolve order parameter gradients across narrow zones where crack surfaces exist. Dynamic remeshing strategies have been invoked in some cases to deal with this issue [137]. Problems associated with proper modeling of large crack velocities via phase field models in the setting of elastodynamics have also been reported [138]. Although the preceding discussion deals with physics of fracture, separation of material in conjunction with shear localization can be addressed via similar principles. A combined phase field model for plastic deformation, adiabatic shear localization, and ductile fracture in metals is described in [139]. Shear band thickness may evolve, which can complicate mesoscale mechanics of plasticity at their interfaces, for example a realistic prescription of a regularization length. Insight into the width and structure of shear bands in ductile metals may be found in analysis and numerical solutions reported in [8, 140, 141]. Another recent diffuse interface approach is applied to shear localization or shear fracture in magnesium [142] and boron carbide [143, 144], with
the latter material known to undergo stress-induced amorphization promoted by compression and shear [145]. The novelty of the recent theory presented in [142, 144]—with initial developments first reported in [146]—is use of Finsler differential geometry to formulate the governing kinematic and equilibrium equations. Essentially, the generalized theory of Finsler-geometric continuum mechanics permits the metric tensor and corresponding local volume elements to depend on the order parameter(s) comprising the state vector of pseudo-Finsler space, enabling a natural coupling between dilatation or volume collapse and shearing modes, without introduction of spurious fitting equations or calibrated parameters. Though not yet undertaken, the same approach could be used to address dilatation in the vicinity of stacking faults or twin boundaries in crystals [147–149]. The first analytical solutions for Finslergeometric continuum theory are derived in [146]; subsequent solutions are also validated versus experimental observations for the above noted crystalline materials in [142, 144, 150, 151]. An alternative Finsler-based thermomechanical theory was applied in numerical simulations of shear localization in metals in [152, 153].
Discussion Capabilities and notable advantages and disadvantages of modeling techniques—sharp interface, continuum defect, and diffuse interface representations—are summarized next. The author’s viewpoint is that all such models have been suitably validated versus experiment or atomic simulation for their intended applications, as has been noted already with supporting references in corresponding sections of the main text of this article and references cited therein. First consider sharp interface models. Sharp interface representations of grain boundaries (GBs) permit explicit modeling of strain and stress concentrations of field variables near mismatched interfaces. No additional constitutive model parameters are required beyond those for the bulk material. The approach is physically realistic for thermoelasticity since properties and hardness tend to vary little from points near the GB to points far in the grain interior. The approach is also considered physically accurate in the sense that atomically sharp boundaries are
J Mater Sci
modeled by jumps in properties, with no artificial smoothing. Similar statements apply for sharp interface models of twin boundaries (TBs). No spurious material parameters generally need to be calibrated to address problems in elastostatics. To address interface motion, however, a kinetic law for GB or TB dynamics must supplement the kinematic description, and a (numerical) scheme must be invoked to track the position of interface(s) as deformation proceeds in time. Implementation of explicit fronttracking often proves challenging. Cohesive failure models are categorized as sharp interface models for material separation. Distinct crack morphologies are suitably captured, as are interactions among cracks and other heterogeneities. Advantages include relatively few parameters for basic crack separation laws and distinct behavior of interfaces and bulk material permitting use of traditional solid continuum finite elements for the bulk. Variations in initial failure properties, which may be extreme in brittle solids, can unfortunately increase the number of parameters needed for realistic predictions. Numerical implementation is usually straightforward. A disadvantage is that fracture paths are often constrained to follow element boundaries. Cohesive finite element sizes must be small enough to resolve behavior over the cohesive length, an issue that may necessitate very fine meshes at correspondingly high computational expense. Next consider continuum defect descriptions of interfaces in crystalline solids. Crystal plasticity can be augmented to explicitly account for defect (e.g., dislocation and disclination) densities in the vicinity of grain and subgrain boundaries. Benefits include physical flexibility, with size effects and microstructure evolution at different resolutions naturally addressed. Regularization associated with higherorder gradient terms may enable mesh-size independence of numerical solutions. Additional parameters must be prescribed and often calibrated rather than determined from first principles. Furthermore, increased model sophistication contributes computational expense and limits availability of analytical solutions for model verification. Pseudo-slip-type twinning models enable reasonably accurate predictions of the onset and evolution of bulk twinning behavior and texture changes. These can usually be implemented in crystal plasticity frameworks with modest effort. Kinetic equations and parameters must be assigned, and effects of twinning on slip must be
included in the flow rule. Calibration and validation of such features is often problematic. Shapes of each twin variant within a local material volume element are not predicted. Continuum damage models can reasonably predict onset and evolution of bulk damage behavior and changes in stiffness. These models can often be implemented in existing simulation frameworks with modest overhead, and no special interfacial elements or node duplications are needed. Kinetic equations and parameters must be assigned for damage and its effects on plastic flow. Calibration is often difficult, and unique specification of all parameters cannot usually be ensured. Complexity increases as multiple scalar damage variables or vector- or tensor-valued damage variables are employed for anisotropic media. Solution fields of discrete cracks or voids and their interactions are not well described. Usually, such models contain no intrinsic regularization length. Finally, consider diffuse interface representations which notably encompass phase field descriptions. Microstructure evolution can often be addressed in a natural way, whereby kinetics of microstructure evolution is motivated by the fundamental principle that a system should seek a minimum energy state. Regularization is obtained as a by-product of gradient terms in the energy functional. Disadvantages include necessary choices of non-unique interpolation functions for property values within interfacial zones and parameters in the energy functional and kinetic laws that may lack obvious physical interpretation. Mesh resolution must be fine enough to resolve field variables and their gradients. Disparate timescales for microstructure evolution and stress dynamics may complicate numerical solution procedures. The diffuse interface approach to modeling twins and twin boundaries enables prediction of detailed, fine-scale twin morphologies for equilibrium states. Diffuse interface models of fracture likewise naturally enable prediction of failure behavior via energy minimization, without recourse to phenomenology usually associated with continuum damage mechanics theories. Relatively, few material parameters or properties are needed. Complex crack morphologies can be predicted by numerical simulations invoking phase field fracture theory. Standard continuum finite elements and their meshes can be used since fields are continuous, but additional degree(s) of freedom are needed to track order parameter(s). A disadvantage is that very fine
J Mater Sci
Table 1 Classes of mesoscopic interfacial models, examples, and general characteristics
Examples (Grain boundary) (Twinning) (Failure) Complexity (usual governing equations) Generality (usual physics addressed) Phenomenology (typical calibration and parameters)
Sharp interface
Continuum defect
Diffuse interface
Discrete GB models (property jumps) Discontinuous deformation gradient Cohesive fracture models Low
Distributed dislocation/disclination
Phase field models of GBs
Pseudo-slip twinning models
Phase field models of twinning
Continuum damage mechanics High
Phase field models of fracture Moderate
Low
High
Moderate
Low
High
Moderate
meshes needed to resolve order parameter gradients across narrow zones (i.e., regularization lengths) where crack surfaces exist. The regularization length is thus often artificially increased in a compromise of physical realism with computational cost. General trends in capabilities are summarized in Table 1, representative of many of the theories developed or implemented in works cited throughout the text. Exceptions to these trends are possible. This article, as implied by its title, has focused on interfacial physics from the perspective of (continuum or mesoscale) mechanics. Applications of sharp and diffuse interface models to broad classes of materials science problems, including those not addressed herein such as solidification, for example, have been discussed in books and other review articles on these subjects [18, 154–156].
Conclusions Models of interfaces in crystalline solids have been categorized, summarized, and evaluated. Emphasis has been given to finite deformation descriptions at the mesoscale, with attention limited to continuum mechanical frameworks, as opposed to molecular statics/dynamics. Particular physical boundaries considered have been grain and subgrain boundaries, twin boundaries, and failure surfaces (e.g., fracture planes). For each type of physical boundary, corresponding models have been grouped into one of three general categories: sharp interface models, continuum defect models, and diffuse interface (e.g., phase field) models. Though moderately lengthy, the
present review of course does not claim to cover all important works, with any notable exclusions unintended. It is hoped that the present article will be a useful aid to researchers beginning studies on the subject and a useful reference for those more experienced, pointing out some key historic and more modern, perhaps overlooked, references wherein further details can be found.
Acknowledgements Much of this paper was written while the author served as a visiting research fellow at Columbia University in the Department of Civil Engineering and Engineering Mechanics of the Fu Foundation School of Engineering and Applied Science in New York, NY, USA. The author acknowledges the courtesy of Dr. WaiChing (Steve) Sun for hosting his sabbatical visit at Columbia University in 2016.
References [1]
[2] [3] [4]
Phillips R (2001) Crystals, defects and microstructures: modeling across scales. Cambridge University Press, Cambridge Rohrer GS (2001) Structure and bonding in crystalline materials. Cambridge University Press, Cambridge Clayton JD (2011a) Nonlinear mechanics of crystals. Springer, Dordrecht Yadav S, Ravichandran G (2003) Penetration resistance of laminated ceramic/polymer structures. Int J Impact Eng 28:557–574
J Mater Sci
[5]
[6] [7]
[8] [9]
[10] [11] [12]
[13]
[14]
[15] [16]
[17]
[18]
[19] [20]
Clayton JD (2015) Modeling and simulation of ballistic penetration of ceramic-polymer-metal layered systems. Math Probl Eng 2015:709498 Brandon DG (1966) The structure of high-angle grain boundaries. Acta Metall 14:1479–1484 Rice JR (1968) Mathematical analysis in the mechanics of fracture. In: Liebowitz H (ed) Fracture: an advanced treatise. Academic Press, New York, pp 191–311 Wright TW (2002) The physics and mathematics of adiabatic shear bands. Cambridge University Press, Cambridge Hughes DA, Hansen N, Bammann DJ (2003) Geometrically necessary boundaries, incidental dislocation boundaries and geometrically necessary dislocations. Scr Mater 48:147–153 Boiko VS, Garber RI, Kosevich AM (1994) Reversible crystal plasticity. AIP Press, New York Christian JW, Mahajan S (1995) Deformation twinning. Prog Mater Sci 39:1–157 Dongare AM, LaMattina B, Irving DL, Rajendran AM, Zikry MA, Brenner DW (2012) An angular-dependent embedded atom method (A-EAM) interatomic potential to model thermodynamic and mechanical behavior of Al/Si composite materials. Modelling Simul Mater Sci Eng 20:035007 Zhigilei LV, Volkov AN, Dongare AM (2012) Computational study of nanomaterials: from large-scale atomistic simulations to mesoscopic modeling. In: Bhushan B (ed) Encyclopedia of nanotechnology. Springer, Berlin, pp 470–480 Abraham FF, Broughton JQ, Bernstein N, Kaxiras E (1998) Spanning the continuum to quantum length scales in a dynamic simulation of brittle fracture. Europhys Lett 44:783–787 Knap J, Ortiz M (2001) An analysis of the quasicontinuum method. J Mech Phys Solids 49:1899–1923 Clayton JD, Chung PW (2006) An atomistic-to-continuum framework for nonlinear crystal mechanics based on asymptotic homogenization. J Mech Phys Solids 54:1604–1639 Chung PW, Clayton JD (2007) Multiscale modeling of point and line defects in cubic crystals. Int J Multiscale Comput Eng 5:203–226 Emmerich H (2003) The diffuse interface approach in materials science: thermodynamic concepts and applications of phase-field models. Springer, Berlin Schoenfeld SE, Wright TW (2003) A failure criterion based on material instability. Int J Solids Struct 40:3021–3037 Wallace DC (2003) Statistical physics of crystals and liquids: a guide to highly accurate equations of state. World Scientific, Singapore
[21]
[22] [23]
[24]
[25] [26]
[27]
[28]
[29]
[30]
[31] [32]
[33] [34] [35]
[36]
[37]
Clayton JD, Tonge A (2015) A nonlinear anisotropic elastic–inelastic constitutive model for polycrystalline ceramics and minerals with application to boron carbide. Int J Solids Struct 64–65:191–207 Bunge H-J (1982) Texture analysis in materials science: mathematical methods. Butterworths, London Randle V, Engler O (2000) Introduction to texture analysis: macrotexture, microtexture, and orientation mapping. Gordon and Breach, Amsterdam Grimmer H, Bollmann W, Warrington DH (1974) Coincident-site lattices and complete pattern-shift lattices in cubic crystals. Acta Crystallogr A 30:197–207 Watanabe T (1984) An approach to grain boundary design for strong and ductile polycrystals. Res Mech 11:47–84 Roters F, Eisenlohr P, Hantcherli L, Tjahjanto DD, Bieler TR, Raabe D (2010) Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: theory, experiments, applications. Acta Mater 58:1152–1211 Clayton JD, Kraft RH, Leavy RB (2012) Mesoscale modeling of nonlinear elasticity and fracture in ceramic polycrystals under dynamic shear and compression. Int J Solids Struct 49:2686–2702 Clayton JD (2013a) Mesoscale modeling of dynamic compression of boron carbide polycrystals. Mech Res Commun 49:57–64 Clayton JD, McDowell DL (2003a) A multiscale multiplicative decomposition for elastoplasticity of polycrystals. Int J Plast 19:1401–1444 Clayton JD, Knap J (2015a) Phase field modeling of directional fracture in anisotropic polycrystals. Comput Mater Sci 98:158–169 Clayton JD (2013b) Nonlinear Eulerian thermoelasticity for anisotropic crystals. J Mech Phys Solids 61:1983–2014 Clayton JD (2014a) Analysis of shock compression of strong single crystals with logarithmic thermoelastic–plastic theory. Int J Eng Sci 79:1–20 Thurston RN (1974) Waves in solids. In: Truesdell C (ed) Handbuch der Physik, vol 4. Springer, Berlin, pp 109–308 Teodosiu C (1982) Elastic models of crystal defects. Springer, Berlin Clayton JD (2014b) Shock compression of metal crystals: a comparison of Eulerian and Lagrangian elastic–plastic theories. Int J Appl Mech 6:1450048 Clayton JD (2015b) Crystal thermoelasticity at extreme loading rates and pressures: analysis of higher-order energy potentials. Mech Lett 3:113–122 Meyers MA, Ashworth E (1982) A model for the effect of grain size on the yield stress of metals. Philos Mag A 46:737–759
J Mater Sci
[38]
[39]
[40]
[41]
[42]
[43]
[44]
[45]
[46]
[47]
[48]
[49]
[50] [51]
[52]
Clayton JD, Schroeter BM, Graham S, McDowell DL (2002) Distributions of stretch and rotation in OFHC Cu. J Eng Mater Technol 124:302–313 Harren SV, Deve HE, Asaro RJ (1988) Shear band formation in plane strain compression. Acta Metall 36:2435–2480 Harren SV, Asaro RJ (1989) Nonuniform deformations in polycrystals and aspects of the validity of the Taylor model. J Mech Phys Solids 37(2):191–232 Clayton JD (2009a) Modeling effects of crystalline microstructure, energy storage mechanisms, and residual volume changes on penetration resistance of precipitatehardened aluminum alloys. Compos B Eng 40:443–450 Needleman A, Asaro RJ, Lemonds J, Peirce D (1985) Finite element analysis of crystalline solids. Comput Methods Appl Mech Eng 52:689–708 Zikry MA, Kao M (1996) Inelastic microstructural failure mechanisms in crystalline materials with high angle grain boundaries. J Mech Phys Solids 44:1765–1798 Ortiz M, Suresh S (1993) Statistical properties of residual stresses and intergranular fracture in ceramic materials. J Appl Mech 60:77–84 Espinos HD, Zavattieri PD (2003a) A grain level model for the study of failure initiation and evolution in polycrystalline brittle materials. Part I theory and numerical implementation. Mech Mater 35:333–364 Pathak S, Michler J, Wasmer K, Kalidindi SR (2012) Studying grain boundary regions in polycrystalline materials using spherical nano-indentation and orientation imaging microscopy. J Mater Sci 47:815–823. doi: 10.1007/s10853011-5859-z Grinfeld M (1991) Thermodynamic methods in the theory of heterogeneous systems. Longman Scientific and Technical, Sussex Abeyaratne R, Knowles JK (1990) On the driving traction acting on a surface of strain discontinuity in a continuum. J Mech Phys Solids 38:345–360 Abeyaratne R, Knowles JK (1991) Kinetic relations and the propagation of phase boundaries in solids. Arch Ration Mech Anal 114:119–154 Read WT, Shockley W (1950) Dislocation models of crystal grain boundaries. Phys Rev 78:275–289 Mughrabi H (1983) Dislocation wall and cell structures and long-range internal stresses in deformed metal crystals. Acta Metall 31:1367–1379 Rezvanian O, Zikry MA, Rajendran AM (2007) Statistically stored, geometrically necessary and grain boundary dislocation densities: microstructural representation and modelling. Proc R Soc Lond A 463:2833–2853
[53]
[54]
[55] [56] [57] [58]
[59]
[60]
[61] [62]
[63]
[64]
[65]
[66]
[67]
[68]
Clayton JD, McDowell DL, Bammann DJ (2006) Modeling dislocations and disclinations with finite micropolar elastoplasticity. Int J Plast 22:210–256 Clayton JD, Bammann DJ, McDowell DL (2004a) Anholonomic configuration spaces and metric tensors in finite strain elastoplasticity. Int J Non Linear Mech 39:1039–1049 Clayton JD (2012a) On anholonomic deformation, geometry, and differentiation. Math Mech Solids 17:702–735 Clayton JD (2014c) Differential geometry and kinematics of continua. World Scientific, Singapore Steinmann P (2015) Geometrical foundations of continuum mechanics. Springer, Berlin Regueiro RA, Bammann DJ, Marin EB, Garikipati K (2002) A nonlocal phenomenological anisotropic finite deformation plasticity model accounting for dislocation defects. J Eng Mater Technol 124:380–387 Clayton JD, McDowell DL, Bammann DJ (2004b) A multiscale gradient theory for elastoviscoplasticity of single crystals. Int J Eng Sci 42:427–457 Admal NC, Po G, Marian J (2017) Diffuse-interface polycrystal plasticity: expressing grain boundaries as geometrically necessary dislocations. Mater Theory 1:1–16 Li JCM (1972) Disclination model of high angle grain boundaries. Surf Sci 31:12–26 Clayton JD, Bammann DJ, McDowell DL (2005) A geometric framework for the kinematics of crystals with defects. Philos Mag 85:3983–4010 Steinmann P (2013) On the roots of continuum mechanics in differential geometry. In: Altenbach H, Eremeyev VA (eds) Generalized continua-from the theory to engineering applications. Springer, Udine, pp 1–64 Clayton JD (2015c) Defects in nonlinear elastic crystals: differential geometry, finite kinematics, and second-order analytical solutions. Z Angew Math Mech ZAMM) 95:476–510 Upadhyay M, Capolungo L, Taupin V, Fressengeas C (2011) Grain boundary and triple junction energies in crystalline media: a disclination based approach. Int J Solids Struct 48:3176–3193 Sun XY, Cordier P, Taupin V, Fressengeas C, Jahn S (2016) Continuous description of a grain boundary in forsterite from atomic scale simulations: the role of disclinations. Philo Mag 96:1757–1772 Gerken JM, Dawson PR (2008) A crystal plasticity model that incorporates stresses and strains due to slip gradients. J Mech Phys Solids 56:1651–1672 Luscher DJ, Mayeur JR, Mourad HM, Hunter A, Kenamond MA (2016) Coupling continuum dislocation
J Mater Sci
[69]
[70]
[71]
[72] [73]
[74]
[75]
[76]
[77]
[78]
[79]
[80] [81]
[82]
[83]
transport with crystal plasticity for application to shock loading conditions. Int J Plast 76:111–129 Clayton JD, Hartley CS, McDowell DL (2014) The missing term in the decomposition of finite deformation. Int J Plast 52:51–76 Clayton JD (2014d) An alternative three-term decomposition for single crystal deformation motivated by non-linear elastic dislocation solutions. Q J Mech Appl Math 67:127–158 Clayton JD, Bammann DJ (2009) Finite deformations and internal forces in elastic–plastic crystals: interpretations from nonlinear elasticity and anharmonic lattice statics. J Eng Mater Technol 131:041201 Toupin RA, Rivlin RS (1960) Dimensional changes in crystals caused by dislocations. J Math Phys 1:8–15 Clayton JD (2009b) A continuum description of nonlinear elasticity, slip and twinning, with application to sapphire. Proc R Soc Lond A 465:307–334 Clayton JD (2009c) A non-linear model for elastic dielectric crystals with mobile vacancies. Int J Non Linear Mech 44:675–688 Abdollahi A, Arias I (2012) Numerical simulation of intergranular and transgranular crack propagation in ferroelectric polycrystals. Int J Fract 174:3–15 Fan D, Chen L-Q (1997) Computer simulation of grain growth using a continuum field model. Acta Mater 45:611–622 Allen SM, Cahn JW (1979) A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening. Acta Metall 27:1085–1095 Abrivard G, Busso EP, Forest S, Appolaire B (2012a) Phase field modelling of grain boundary motion driven by curvature and stored energy gradients. Part I: theory and numerical implementation. Philos Mag 92:3618–3642 Abrivard G, Busso EP, Forest B, Appolaire S (2012b) Phase field modelling of grain boundary motion driven by curvature and stored energy gradients. Part II: application to recrystallisation. Philos Mag 92:3643–3664 James RD (1981) Finite deformation by mechanical twinning. Arch Ration Mech Anal 77:143–176 Zanzotto G (1996) The Cauchy–Born hypothesis, nonlinear elasticity and mechanical twinning in crystals. Acta Crystallogr A 52:839–849 Bilby BA, Crocker AG (1965) The theory of the crystallography of deformation twinning. Proc R Soc Lond A 288:240–255 Barton NR, Winter NW, Reaugh JE (2009) Defect evolution and pore collapse in crystalline energetic materials. Modelling Simul Mater Sci Eng 17:035003
[84]
[85]
[86]
[87]
[88]
[89]
[90]
[91]
[92]
[93]
[94]
[95]
[96]
[97]
[98]
[99]
Clayton JD, Becker R (2012) Elastic–plastic behavior of cyclotrimethylene trinitramine single crystals under spherical indentation: modeling and simulation. J Appl Phys 111:063512 Bhattacharya K (2003) Microstructure of martensite: why it forms and how it gives rise to the shape-memory effect. Oxford University Press, New York Hou TY, Rosakis P, LeFloch P (1999) A level-set approach to the computation of twinning and phase-transition dynamics. J Comput Phys 150:302–331 Chin GY, Hosford WF, Mendorf DR (1969) Accommodation of constrained deformation in FCC metals by slip and twinning. Proc R Soc Lond A 309:433–456 Van Houtte P (1978) Simulation of the rolling and shear texture of brass by the Taylor theory adapted for mechanical twinning. Acta Metall 26:591–604 Kalidindi SR (1998) Incorporation of deformation twinning in crystal plasticity models. J Mech Phys Solids 46:267–290 Staroselsky A, Anand L (1998) Inelastic deformation of polycrystalline face centered cubic materials by slip and twinning. J Mech Phys Solids 46:671–696 Clayton JD (2010a) Modeling finite deformations in trigonal ceramic crystals with lattice defects. Int J Plast 26:1357–1386 Mirkhani H, Joshi SP (2014) Mechanism-based crystal plasticity modeling of twin boundary migration in nanotwinned face-centered-cubic metals. J Mech Phys Solids 68:107–133 Clayton JD, Knap J (2011a) A phase field model of deformation twinning: nonlinear theory and numerical simulations. Physica D 240:841–858 Hu SY, Henager CH, Chen L-Q (2010) Simulations of stress-induced twinning and de-twinning: a phase field model. Acta Mater 58:6554–6564 Heo TW, Wang Y, Bhattacharya S, Sun X, Hu S, Chen L-Q (2011) A phase-field model for deformation twinning. Philos Mag Lett 91:110–121 Bhattacharya K (1993) Comparison of the geometrically nonlinear and linear theories of martensitic transformation. Contin Mech Thermodyn 5:205–242 Clayton JD, Knap J (2011b) Phase field modeling of twinning in indentation of transparent single crystals. Modelling Simul Mater Sci Eng 19:085005 Clayton JD, Knap J (2013) Phase field analysis of fracture induced twinning in single crystals. Acta Mater 61:5341–5353 Hildebrand FE, Miehe C (2012) A phase field model for the formation and evolution of martensitic laminate microstructure at finite strains. Philos Mag 92:4250–4290
J Mater Sci
[100] Clayton JD, Knap J (2015b) Nonlinear phase field theory for fracture and twinning with analysis of simple shear. Philos Mag 95:2661–2696 [101] Kondo R, Tadano Y, Shizawa K (2014) A phase-field model of twinning and detwinning coupled with dislocation-based crystal plasticity for HCP metals. Comput Mater Sci 95:672–683 [102] Li X, Wei Y, Lu L, Lu K, Gao H (2010) Dislocation nucleation governed softening and maximum strength in nano-twinned metals. Nature 464:877–880 [103] Yang L, Dayal K (2010) Formulation of phase-field energies for microstructure in complex crystal structures. Appl Phys Lett 96:081916 [104] Agrawal V, Dayal K (2015a) A dynamic phase-field model for structural transformations and twinning: regularized interfaces with transparent prescription of complex kinetics and nucleation. part I: formulation and one-dimensional characterization. J Mech Phys Solids 85:270–290 [105] Agrawal V, Dayal K (2015b) A dynamic phase-field model for structural transformations and twinning: regularized interfaces with transparent prescription of complex kinetics and nucleation. part II: two-dimensional characterization and boundary kinetics. J Mech Phys Solids 85:291–307 [106] Schultz MC, Jensen RA, Bradt RC (1994) Single crystal cleavage of brittle materials. Int J Fract 65:291–312 [107] Clayton JD (2008) A model for deformation and fragmentation in crushable brittle solids. Int J Impact Eng 35:269–289 [108] Antoun T (2003) Spall fracture. Springer, New York [109] Dugdale DS (1960) Yielding of steel sheets containing slits. J Mech Phys Solids 8:100–104 [110] Barenblatt GI (1962) The mathematical theory of equilibrium cracks in brittle fracture. Adv Appl Mech 7:55–129 [111] Needleman A (1987) A continuum model for void nucleation by inclusion debonding. J Appl Mech 54:525–531 [112] Xu X-P, Needleman A (1993) Void nucleation by inclusion debonding in a crystal matrix. Modelling Simul Mater Sci Eng 1:111–132 [113] Clayton JD, McDowell DL (2004) Homogenized finite elastoplasticity and damage: theory and computations. Mech Mater 36:799–824 [114] Xu X-P, Needleman A (1994) Numerical simulations of fast crack growth in brittle solids. J Mech Phys Solids 42:1397–1434 [115] Clayton JD (2005a) Dynamic plasticity and fracture in high density polycrystals: constitutive modeling and numerical simulation. J Mech Phys Solids 53:261–301 [116] Clayton JD (2005b) Modeling dynamic plasticity and spall fracture in high density polycrystalline alloys. Int J Solids Struct 42:4613–4640
[117] Vogler TJ, Clayton JD (2008) Heterogeneous deformation and spall of an extruded tungsten alloy: plate impact experiments and crystal plasticity modeling. J Mech Phys Solids 56:297–335 [118] Espinosa HD, Zavattieri PD (2003b) A grain level model for the study of failure initiation and evolution in polycrystalline brittle materials. Part II: numerical examples. Mech Mater 35:365–394 [119] Kraft RH, Molinari JF (2008) A statistical investigation of the effects of grain boundary properties on transgranular fracture. Acta Mater 56:4739–4749 [120] Foulk JW, Vogler TJ (2010) A grain-scale study of spall in brittle materials. Int J Fract 163:225–242 [121] Krajcinovic D (1996) Damage mechanics. Elsevier, Amsterdam [122] Voyiadjis GZ, Kattan PI (2005) Damage mechanics. CRC Press, Boca Raton [123] Kachanov LM (1958) Time of the rupture process under creep conditions. Isv Akad Nauk SSR Otd Tekh Nauk 8:26–31 [124] Bammann DJ, Solanki KN (2010) On kinematic, thermodynamic, and kinetic coupling of a damage theory for polycrystalline material. Int J Plast 26:775–793 [125] Del Piero G, Owen DR (1993) Structured deformations of continua. Arch Ration Mech Anal 124:99–155 [126] Clayton JD (2006) Continuum multiscale modeling of finite deformation plasticity and anisotropic damage in polycrystals. Theor Appl Fract Mech 45:163–185 [127] Clayton JD, McDowell DL (2003b) Finite polycrystalline elastoplasticity and damage: multiscale kinematics. Int Solids Struct 40:5669–5688 [128] Clayton JD (2010b) Deformation, fracture, and fragmentation in brittle geologic solids. Int J Fract 163:151–172 [129] Bammann DJ, Aifantis EC (1989) A damage model for ductile metals. Nucl Eng Des 116:355–362 [130] Aslan O, Cordero NM, Gaubert A, Forest S (2011) Micromorphic approach to single crystal plasticity and damage. Int J Eng Sci 49:1311–1325 [131] Jin YM, Wang YU, Khachaturyan AG (2001) Three-dimensional phase field microelasticity theory and modeling of multiple cracks and voids. Appl Phys Lett 79:3071–3073 [132] Karma A, Kessler DA, Levine H (2001) Phase-field model of mode III dynamic fracture. Phys Rev Lett 87:045501 [133] Del Piero G, Lancioni G, March R (2007) A variational model for fracture mechanics: numerical experiments. J Mech Phys Solids 55:2513–2537 [134] Clayton JD, Knap J (2014) A geometrically nonlinear phase field theory of brittle fracture. Int J Fract 189:139–148
J Mater Sci
[135] Clayton JD, Knap J (2016) Phase field modeling of coupled fracture and twinning in single crystals and polycrystals. Comput Methods Appl Mech Eng 312:447–467 [136] Spatschek R, Brener E, Karma A (2011) Phase field modeling of crack propagation. Philos Mag 91:75–95 [137] Borden MJ, Verhoosel CV, Scott MA, Hughes TJR, Landis CM (2012) A phase-field description of dynamic brittle fracture. Comput Methods Appl Mech Eng 217:77–95 [138] Agrawal V, Dayal K (2017) Dependence of equilibrium Griffith surface energy on crack speed in phase-field models for fracture coupled to elastodynamics. Int J Fract 207:243–249 [139] McAuliffe C, Waisman H (2015) A unified model for metal failure capturing shear banding and fracture. Int J Plast 65:131–151 [140] Wright TW, Ockendon H (1992) A model for fully formed shear bands. J Mech Phys Solids 40:1217–1226 [141] Wright TW, Walter JW (1996) The asymptotic structure of an adiabatic shear band in antiplane motion. J Mech Phys Solids 44:77–97 [142] Clayton JD (2017a) Finsler geometry of nonlinear elastic solids with internal structure. J Geom Phys 112:118–146 [143] Clayton JD (2014e) Phase field theory and analysis of pressure-shear induced amorphization and failure in boron carbide ceramic. AIMS Mater Sci 1:143–158 [144] Clayton JD (2017b) Generalized finsler geometric continuum physics with applications in fracture and phase transformations. Z Angew Math Phys (ZAMP) 68:9 [145] Clayton JD (2012b) Towards a nonlinear elastic representation of finite compression and instability of boron carbide ceramic. Philos Mag 92:2860–2893
[146] Clayton JD (2016a) Finsler-geometric continuum mechanics. Technical Report ARL-TR-7694, US Army Research Laboratory, Aberdeen Proving Ground MD [147] Kenway PR (1993) Calculated stacking-fault energies in aAl2 O3 . Philos Mag B 68:171–183 [148] Clayton JD (2010c) Modeling nonlinear electromechanical behavior of shocked silicon carbide. J Appl Phys 107:013520 [149] Clayton JD (2011b) A nonlinear thermomechanical model of spinel ceramics applied to aluminum oxynitride (AlON). J Appl Mech 78:011013 [150] Clayton JD (2016b) Finsler-geometric continuum mechanics and the micromechanics of fracture in crystals. J Micromech Mol Phys 1:164003 [151] Clayton JD (2017) Finsler-geometric continuum dynamics and shock compression. Int J Fract. doi:10.1007/s10704017-0211-5 [152] Saczuk J (1996) Finslerian foundations of solid mechanics. Polskiej Akademii Nauk, Gdansk [153] Stumpf H, Saczuk J (2000) A generalized model of oriented continuum with defects. Z Angew Math Mech (ZAMM) 80:147–169 [154] Sethian JA (1999) Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science. Cambridge University Press, Cambridge [155] Steinbach I (2009) Phase-field models in materials science. Modelling Simul Mater Sci Eng 17:073001 [156] Chen L-Q (2002) Phase-field models for microstructure evolution. Ann Rev Mater Res 32:113–140