A Hybrid Multi-Scale Model of Crystal Plasticity for Handling Stress Concentrations

Microstructural effects become important at regions of stress concentrators such as notches, cracks and contact surfaces. A multiscale model is presented that efficiently captures microstructural details at such critical regions. The approach is based on a multiresolution mesh that includes an explicit microstructure representation at critical regions where stresses are localized. At regions farther away from the stress concentration, a reduced order model that statistically captures the effect of the microstructure is employed. The statistical model is based on a finite element representation of the orientation distribution function (ODF). As an illustrative example, we have applied the multiscaling method to compute the stress intensity factor K I around the crack tip in a wedge-opening load specimen. The approach is verified with an analytical solution within linear elasticity approximation and is then extended to allow modeling of microstructural effects on crack tip plasticity.


Introduction
Efficient micro-scale modeling tools are needed to compute microstructure-dependent properties of advanced structural alloys used in aerospace, naval and automotive applications [1].Microstructural effects become an important consideration in regions of stress concentrations such as notches, cracks and contact surfaces.While the crystal plasticity finite element (CPFE) method has emerged as an effective tool for simulating the mechanical response of aggregates of a few hundred metallic crystals, the simulation of 'macro-scale' components that contain millions of grains is a challenging task even when using current state-of-the-art supercomputers.Multiscale methods developed to address this problem can be broadly classified into two categories: concurrent and hierarchical methods.In hierarchical methods, the objective is to obtain material properties of a coarse-scale model by simulating a fine-scale model that contains features smaller than those allowed by the coarse-scale mesh resolution.A popular example is the computational homogenization approach [2,3] where volume averaging schemes are employed to extract macroscale properties.
In concurrent methods, both fine-and coarse-scale problems are solved together in a single domain using multi-resolution meshes.Most popular among these methods are quasicontinuum methods, where the resolution of meshes from atomistic to continuum scales has been considered to study crack propagation [4].In QC methods, the treatment of interfaces between fine-and coarse-scale regions becomes problematic when multi-physical domains (e.g., crystal plasticity at the fine scale and continuum plasticity at the coarse scale) are considered.For example, Dunne et al. (2007) [5] employed a combined continuum crystal plasticity finite element approach for modeling a three-point bend test.One of the potential difficulties identified in the model was the inconsistency of yield stress at the interface between the domains.To minimize the effect, the choice of the critical resolved shear stress in the crystal plasticity model was made such that polycrystal yielding occurred at a value close to the continuum yield stress.
In this paper, the interface problem is solved through the use of a single physical model at both the fine and coarse scales.This is achieved (i) by explicitly resolving single-crystal deformation at the fine scale (crack tip) using a crystal plasticity model and (ii) by employing a hierarchical multiscale model at the coarse scale (far region), where polycrystal aggregates are homogenized using the same crystal plasticity model.While the method combines the best aspects of both hierarchical and concurrent methods, it is still computationally expensive to model polycrystalline aggregates at every coarse-scale material point.
An alternate class of schemes has been developed in recent years that allows the representation of microstructure using probabilistic descriptors.The simplest of these descriptors is the one-point probability measure, the orientation distribution function (A(g)), which quantifies the volume fractions of crystals in the orientation space (g).Under an applied deformation, texturing is simulated by numerically evolving the orientation distribution function (ODF) (rather than the microstructure itself) using conservation laws [6].A variety of ODFs can be represented over a finite element (FE) mesh of the orientation space that has significantly lower degrees of freedom when compared to a microstructure FE mesh [7][8][9].For example, individual crystals close to the crack tip can be modeled using an ODF that is a delta function in the orientation space.Farther from the notch tip, smoothly varying ODFs can be modeled over an uniform grid.The overall multi-scale model is shown in Figure 1.Here, the transition from the microscale to the macroscale is managed through adaptive meshing methods where element lengths can vary from micrometers at the critical regions to as large as several centimeters at the scale of the macroscale component.As an illustrative example, we have applied the multiscaling method to a well-studied problem in fracture mechanics: the computation of the stress intensity factor K I around the crack tip in a wedge-opening load specimen [10].The multiscale model is subsequently used to analyze the effect of the microstructure at the crack tip when crystal plasticity effects are included.The paper is arranged as follows: In Section 2, the statistical representation of the microstructure is briefly explained, followed by the single-crystal constitutive model.In Section 3, the macro-scale finite element solution procedure is explained.In Section 4, an illustrative example of the approach is provided, followed by conclusions in Section 5.

Representation of the ODF
The orientation distribution function (ODF, A(g)) gives the probability density of finding an orientation g in the microstructure.The ODF satisfies the following conservation equations at all times during deformation: where dg is a differential volume element (the invariant measure) of the orientation space.In addition to the above constraints, the orientation space corresponding to all possible g's must satisfy the crystallographic symmetries of the chosen system (FCC, HCP, etc.).The complete orientation space of a polycrystal can be reduced to a smaller subset, called the fundamental region, as a consequence of crystal symmetries.Within the fundamental region, each crystal orientation is represented uniquely by the coordinate g, the parametrization for the rotation (e.g., Euler angles, Rodrigues vector etc.).The ODF (A(g)) can be represented as a probability density function over the fundamental region of orientation space.In this work, a finite element (FE) mesh is used to model the fundamental region, and the ODF is defined at the nodal points of this mesh [8,9].The probability values between nodal points are obtained as a result of interpolation using finite element shape functions.

Probability Update in Finite Element Spaces
The probabilities are evolved from time t = 0 from an initial ODF that satisfies the conservation equations in Equation (1) using a Lagrangian finite element approach [6,7].The initial orientation g o of a crystal reorients during deformation and maps to a new orientation g t at time t.Simultaneously, the finite element mesh of fundamental region M g deforms with nodes located at g o moving to new locations g t .We assume that the mapping from g o to g t is invertible.
The ODF A(g t ) represents the probability density of crystals with orientation g t at time t.The evolution of the ODF is given by the conservation Equation (1) as: where .Using the Jacobian, a map of the current mesh (at time t) to the reference mesh (at t = 0) can be made: The quantity written as Â(g o , t) is the volume density A(g t ) plotted over the corresponding orientation (g o ) in the initial mesh.Thus, Â(g o , t) gives the Lagrangian representation of the current ODF in the initial mesh M g o .If the integrand is continuous, a localized relationship of the following form can be used to update the ODF at any time t:

ODF for Planar Polycrystals
For simplicity, consider planar polycrystals characterized by a two-dimensional rotation R that relates the local crystal lattice frame to the reference sample frame.A parametrization of the associated rotation group is, R = Icos(g) − Esin(g) (5) where g is the angle between the crystal and sample axes, E is the two-dimensional alternator (E 11 = E 22 = 0, E 12 = −E 21 = 1) and I is the identity tensor.A general planar crystal with symmetry under rotations through π is considered here.Under the symmetry, crystal orientations can be described uniquely by parameters drawn from a simply connected fundamental region [a, a + π).
Out of convenience, we will restrict the choice of fundamental regions to the interval closest to the origin (−π/2, π/2).Due to symmetry, the orientation π/2 is exactly the same as orientation −π/2.
The symmetry constraint on the ODF is enforced by using periodic boundary conditions in the finite element mesh.Figure 2 gives an idea of how the approach works for a one-dimensional fundamental region that is represented using two-noded finite elements with linear interpolation.Here, the Jacobian is simply the ratio of element lengths, i.e., current length divided by the initial length.If the element length decreases over time, the probability density has to increase based on Equation ( 4) to maintain normalization of the ODF.Note that the integrand in Equation (3) needs to be continuous for the localization relationship to be valid.Thus, it is implied that J F (g o , t) needs to be continuous.For computing J F (g o , t), ∂t is calculated from the elasto-plastic constitutive model (described next).At time t during deformation, the new positions (g t ) of nodes in the fundamental region mesh are computed using the reorientation velocities v from the constitutive model.The expression for v is obtained by taking a derivative of relation Equation ( 5): where Ω is the spin tensor defined as Ω = Ṙe R e T .Here, R e is evaluated through the polar decomposition of the elastic part of deformation gradient (F e = R e U e ) as computed from the constitutive model (described in Section 2.3).The reorientation velocity v = ∂g t ∂t is computed at each nodal point in the mesh and the change in orientation ∆g = g t − g o is then calculated and stored at the nodal points in the fundamental region.The Jacobian (J F (g o , t)) is then computed using finite element shape functions.
The average stress for the microstructure is obtained by averaging the single crystal stresses σ over the ODF as follows: where σ(g o , t) is the stress (σ(g t ) = σ(g o + ∆g)) plotted over the corresponding orientation (g o ) in the initial ODF mesh.4), which ensures that the normalization constraint (Equation ( 2)) is met in the reoriented mesh.

Constitutive Modeling
We employ the classical single-crystal plasticity theory [11] based on the notion that plastic flow takes place through slip on prescribed slip systems.For a material with α = 1, . . ., N slip systems defined by ortho-normal vector pairs (m α 0 , n α 0 ) denoting the slip direction and slip plane normal respectively at time t = 0, the constitutive equations relate the following basic fields (all quantities expressed in crystal lattice coordinate frame): the deformation gradient F defined with respect to the initial undeformed crystal, which can be decomposed into elastic and plastic parts as F = F e F p (with det (F p ) = 1), the Cauchy stress σ and the slip resistances s α > 0. In the constitutive equations to be defined below, the Green elastic strain measure E e = 1  2 F e T F e − I is defined on the intermediate configuration (plastically deformed, unstressed configuration) utilized.The conjugate stress measure is then defined as T = detF e (F e ) −1 σ(F e ) −T .
The constitutive relation, for stress, is given by T = L e [E e ] where L e is the fourth-order anisotropic elasticity tensor.It is assumed that deformation takes place through dislocation glide, and the evolution of the plastic velocity gradient is given by: where S α 0 = m α 0 ⊗ n α 0 is the Schmid tensor and γα is the plastic shearing rate on the α-th slip system.The resolved stress on the α-th slip system is given by τ α = T • S α 0 .The evolution of slip system resistance is given by the following expression: where, ) a (no sum on β) (10) Here, the slip system hardening term (h αβ ) includes latent hardening through parameter q.In this term, h β o is the hardening coefficient of slip system β, and s β s is the saturation resistance of slip system β.
The resolved shear stress τ α = T • S α 0 is taken to attain a critical value s α on the systems where slip occurs with the plastic shearing rate on the α-th slip system γα > 0. Further, the resolved shear stress does not exceed s α on the inactive systems with γα = 0.A potentially active set PA of slip systems can be identified initially based on the trial resolved stress as the systems with |τ α tr | − s α > 0, where τ α tr = L e [E e tr ] • S α 0 with E e tr = 1 2 (F e tr ) T F e tr − I and F e tr = FF p−1 n , where F p n is the plastic deformation gradient at the previous time step.
During plastic flow, the active systems are assumed to follow the consistency condition: |τ α | = s α .The increment in shearing rates ∆γ β at each time step is obtained by solving the following equation: where, α, β ∈ PA, B β = (S β 0 ) T (F e tr ) T F e tr + (F e tr ) T F e tr S β 0 .A system of equations is obtained of the following form, where, Following this step, the plastic and elastic parts of the deformation gradient are updated for each orientation.The slip resistances are also updated at the end of the time step.

Finite Element Algorithm
The equilibrium equation is expressed in the Lagrangian framework as: where ∇ 0 is the divergence described in the initial reference configuration (t = 0) and f denotes the body forces.The polycrystal Piola-Kirchhoff-I stress, P , is defined as: For any kinematically admissible test function ũ, the weak form of the virtual work equation with applied surface traction λ can be written as: The Newton-Raphson iterative scheme with a line search procedure is employed to solve the weak form while using a formulation that uses a combination of triangle and quadrilateral elements: The linearization process of PK-I stress is given by: δσ can be expressed as: where δF e is obtained as: Here, δT can be obtained as:  The average stress and tangent modulus at the macro-scale level are obtained by averaging over the ODF as given in Equation (7).At the micro-scale, the tangent modulus and stress are exactly computed in the elements within individual crystals.

Numerical Results
A specific crystal geometry with two slip systems at orientation −π/6 and +π/6 is considered [12].In the 2-D elastic stiffness matrix, c 11 = 2 GPa and c 12 = c 44 = 1 GPa.The calculation is done under plane stress conditions.The parameters in the hardening law are taken to be as follows: 4 and a = 2.The single edge notch configuration (also called the wedge-opening load (WOL) specimen) is depicted in Figure 3.The calculation is only performed using one half of the model for simplicity.A 5000 N shear force p is applied on the left edge.Specific dimensions are a/w = 0.5 and h/w = 0.6, as shown in Figure 1.All of the figures are plotted using non-dimensional variables, K I b √ w p and r/w. a is the notch length; w is the width; 2h is the height of the specimen; b is unit thickness; for the plane stress problem, b = 1 m.The solution leads to a stress singularity at the notch tip.There is a wealth of literature on the use of special elements to capture this crack tip singularity, the most popular of which are the quarter point elements that capture the 1/ √ r dependence of the crack tip singularity in linear elasticity or the 1/r dependence in perfect plasticity.However, the singularity of the crack tip solution in crystal plasticity is not known analytically, which makes the design of such elements non-trivial.For simplicity, conventional finite elements (the combination of three-noded triangular and four-noded rectangular elements) are employed in these simulations.Meaningful values of crack tip stress intensity factors can be obtained using this approach with the computational procedures demonstrated in Chan (1970) [10].We only consider the Mode I crack in our paper.The stress intensity factor calculation follows Chan's paper [10].Stress intensity K * I from the displacement-based method: 2 , v * is the computed nodal displacement in the y direction.The K * I is obtained by substituting a calculated nodal point displacement v * at point (r, θ).If the v * were the exact theoretical value, then the value of K * I will also be exact.The above equation is for the plane strain situation; for plane stress cases, just replace ν with ν/(1 + ν) for f (θ, ν).
The stress intensity factor K * I from the stress-based method was also used for comparison: . This theoretical result is obtained from the collocation method described in [13].
The J integral is calculated as: The Γ is an arbitrary curve surrounding the crack tip.W is the strain energy density, while T is the traction vector defined by the out normal along Γ, T i = σ ij n j , and u is the displacement vector.The integral is evaluated counterclockwise around the crack tip.The reported J value is twice the half plane's J value.Strain energy density W is found from PK-I stress and deformation gradient F, W = Ω 0 P : δFdV.The shear force on the left boundary of the model is applied over 20 steps to get the value of p.

Linear Elastic Simulations
The first few tests were performed by setting the plastic shearing strain to zero ∆γ = 0 which leads to a fully-elastic model.In the convergence results shown below, the displacement on the crack surface (θ = π) and stress σ yy on the θ = 0 plane (away from the crack in the positive xdirection) are shown.First, a convergence test is performed on different meshes by progressively refining the element size in the microstructure; see Figure 4.In this simulation, a fully-isotropic model is employed without any microstructure.The test shows good convergence with the increasing number of elements, as shown in Figure 5, and the results are similar to those reported in the finite element model of [10].The analytical result is from [13].
A similar convergence test was performed by altering the number of grains in the vicinity of the crack tip (and 835 elements in each case; see Figure 6).The zero grain case corresponds to the isotropic elastic case (with a fully-random ODF).This corresponds to an isotropic stiffness matrix with c 11 = 2.25 GPa, c 12 = 0.75 GPa and c 44 = 0.75 GPa.The grain convergence study allows one to choose the size of the critical region that needs to be modeled using microstructural grids.From the results shown in Figure 7, there is a trend towards convergence with the increasing number of grains.As the number of grains increases from 42 to 62, little additional accuracy is realized.This implies that a 42-grain case is sufficient to address this particular problem.Beyond this point, the ODF model may be employed.Further, a quick comparison of the zero grain case (isotropic solution) with the 62 grain case result in Figure 7 clearly depicts the role of the microstructure in determining the stress concentration at the crack tip.When r/w is smaller than 6 × 10 −3 , stress concentrations deviate from the isotropic case.While r/w is larger than 6 × 10 −3 , the microstructure does not play a significant role, and the result converges back to the isotropic solution.
Subsequently, three models are compared: (i) a model with isotropic properties; (ii) a model with every element having an underlying ODF (called an 'all-ODF' model) with no explicit microstructure; (iii) a multiscale model with explicitly-modeled crystal orientations (called 'Elastic-Ot'model) near the crack tip and ODF representation at the far field.The detailed comparison between these models is shown in Figure 8. Firstly, the all-ODF model (with the initial random ODF) has the exact same result as the isotropic elastic model as expected.The model with the explicit microstructure has larger variations in stress concentration near the crack tip. Figure 9 additionally depicts the resolved shear stress in the two slip systems for the elastic model as overlaid on the grains showing significant variability in stresses from grain to grain, stressing the need for a microstructure model at critical regions.

Elasto-Plastic Simulations
In subsequent simulations, notch tip plasticity is turned on.All test cases are based on meshes with 835 elements and 62 grains cover the notch tip area.The inclusion of notch tip plasticity leads to a larger displacement at the crack tip (and also smaller stress near the crack tip).This is seen in Figure 10 where the displacements and stresses are directly compared to the elastic case.These results are subsequently shown in terms of the non-dimensional stress intensity factor in Figure 11.Note that on the θ = π plane, the stress intensity factor Equation (25) as computed from the elasticity-based formula will be larger in the plastic case due to larger displacements.On the other hand, the smaller stress σ yy on the θ = 0 plane will lead to a smaller K I value based on Equation (26).These are simply artifacts of using an elasticity-based analytical solution to calculate the stress intensity factor.Due to the larger displacement and smaller stresses around crack tips, elasto-plastic results show altogether different K I trends when compared to the fully-elastic case in Figure 8.
Next, comparing an all-ODF model with plasticity (zero grains) and one with the microstructure explicitly modeled, the differences between these two cases (shown in Figure 11) are much more pronounced than the fully-elastic case (that did a similar comparison between no grains and the 62-grain case; see Figure 7).The inclusion of the microstructure significantly affects the calculated stress intensity factor, which is important during the assessment of safety factors during component design.Figure 12a,b additionally depicts the resolved shear stress in the two slip systems for the elasto-plastic model as overlaid on the grains.This again shows significant variability in stresses from grain to grain, but stresses are smaller than the fully-elastic case shown in Figure 9 as expected.As addressed in McMeeking [14] and Rice and Johnson [15], the crack tip plastic blunting can also be explained by exploiting the crack tip opening stresses.However, instead of using the analytical method, we will employ the J-integral approach next to more reliably compare the elastic and elastoplastic cases.

J Integral Calculation
In order to numerically evaluate the J integral, several loops are drawn around the crack tip as shown in Figure 13.The inner-most curve is within only one grain, while the outmost curve circumscribes the entire microstructure.The results of the J integral calculations are shown in Figure 14.The J integrals of multiscale models (with microstructure) are normalized with respect to that of a fully-homogenized model (with no grains).The x-coordinate is the mean distance values of points on different integral curves.Figure 14a shows the J values of the elastic model.The elastic model with the microstructure has lower J values at the crack tip, but converges to a homogenized model at the far field.Figure 14b shows the J integral of the elastoplastic model.This model also converges to the homogenized model result at the far field; however, the near field J values oscillate because of the microstructural influence.This trend is similar to those shown in Figures 8 and 11 in that the model with the microstructure has more prominent oscillation of stress intensities around the crack tip in the case of an elastoplastic model.
In the elastic case, the ODF does not evolve with deformation.However, in the elastoplastic case, the texture evolves close to the crack tip with the specimen deformation.Texture evolves in the region modeled explicitly with the microstructure, at the same time, and Figure 15 shows the case of ODF evolution in the homogenized region.

Conclusions
The proposed multi-scale model combines the computational efficiency of statistical approaches with the accuracy of CPFEat critical regions.The approach can provide detailed stress and displacement results around critical regions such as crack tips and localized loads (e.g., indentation tests) while also providing good estimates of far field stresses.As an illustrative example, we applied the multiscaling method to compute the stress intensity factor K I around the crack tip in a wedge-opening load specimen [10].The approach was verified with an analytical solution within linear elasticity and has been extended to allow the modeling of crack tip plasticity.We are currently extending this framework to allow modeling of arbitrary crack paths at critical regions within the multiscale approach.

Figure 1 .
Figure 1.The configuration of the multi-scale model.The model is divided into two scales.The microscale is the critical areas containing crystals with orientations.The macroscale element has an orientation distribution function (ODF) in each integration point representing hundreds to thousands of crystals.

Figure 2 .
Figure2.Probability update scheme in finite element (FE) space: During deformation, the nodal points (g) of the FE mesh are moved to reflect the reorientation (∆g) of crystals.The new ODF is obtained using Equation (4), which ensures that the normalization constraint (Equation (2)) is met in the reoriented mesh.

Figure 3 .
Figure 3.The single-edge notch configuration with width w, height 2h and shear force p is applied on the left edge.Only the upper half is modeled.

Figure 5 .
Figure 5. Grid convergence test with 530 elements, 835 elements, 1033 elements and 1225 elements of 62 grains on the elastic all-ODF model.The analytical line is from the collocation method.(a) Non-dimensional K I from the displacement method; (b) non-dimensional K I from the stress method.Both results show convergence at 835 elements.

Figure 6 .Figure 7 .Figure 8 .Figure 9 .
Figure 6.Grids with increasing number of grains near the crack tip (22, 42 and 62 grains) with 835 elements.Different numbers of grains are used to test the critical area that needs to be covered by

Figure 13 .Figure 14 .
Figure13.J integral curves around the crack tip: the red lines are curves; blue lines are crystals; and black ones are FE elements.The innermost curve is only within one crystal.Other curves cover larger areas than the first one.Each edge of the curve is divided into 10 parts.This means that if the curve has five edges, and there are 50 points to be integrated for the J integral.All distances are in mm.

Figure 15 .
Figure 15.The final texture returned by in the plastic model at the macro-scale elements.The final ODFs at different integration points are different from each other; while for the elastic result, the final ODFs at differentintegration points are the same.

where G is the reorientation gradient given as G(g o , t) = ∂g t ∂g o
dg o represents the volume element in the undeformed (initial) ODF mesh (M g o ), which becomes volume element dg t at time t.A Jacobian J F (g o , t) = det(G) gives the ratio of elemental volumes such that dg t = J F (g o , t)dg o ,