A Large-Deformation Gradient Damage Model for Single Crystals Based on Microdamage Theory

: This work aims at the uniﬁcation of the thermodynamically consistent representation of the micromorphic theory and the microdamage approach for the purpose of modeling crack growth and damage regularization in crystalline solids. In contrast to the thermodynamical representation of the microdamage theory, micromorphic contribution to ﬂow resistance is deﬁned in a dual fashion as energetic and dissipative in character, in order to bring certain clarity and consistency to the modeling aspects. The approach is further extended for large deformations and numerically implemented in a commercial ﬁnite element software. Speciﬁc numerical model problems are presented in order to demonstrate the ability of the approach to regularize anisotropic damage ﬁelds for large deformations and eliminate mesh dependency.


Introduction
Highly anisotropic behavior of single crystals necessitates the consideration of microstructurally motivated life time estimation models. Considering the general use of single crystals in complex geometries under harsh environments such as turbine blades, modeling of crack growth due to the coupling of plasticity and damage becomes essential [1,2]. The chosen constitutive model should satisfactorily demonstrate the extremely anisotropic and nonlinear behavior of single crystals under severe thermomechanical loading conditions. There are many modeling attempts in the literature which associate damage localization to crystallographic planes and inelastic deformations [3][4][5].
Conventional rate-independent strain based approaches do not effectively model the localization of deformation due to damage within a finite thickness, since they lack an intrinsic length scale. A theory which includes a material length scale is necessary for modeling such an inherently size-dependent phenomenon. Earliest attempts of modeling introducing a length scale are in the plasticity theory [6,7]. The proposed approach is further extended and numerically implemented in several works [8][9][10] and established as so-called gradient plasticity theories with the remarkable works by Engelen, Peerlings, Geers and Forest [11][12][13][14][15]. Based on principle of virtual power, a thermodynamically consistent formulation considering nonlocal terms energetic, i.e., not dissipative in character, is presented by Gurtin and Anand [16] and a micromorphic version of the theory is also demonstrated in Anand et al. [17]. In the same spirit, there are many regularization attempts in the literature for strain softening behavior due to damage [18]. Such a micromorphic continuum based regularization technique has also been proposed for single crystals in order to eliminate the mesh-dependency problem of the results of such simulations of finite elements by Aslan et al. [19]. However, those attempts exclude the energetic nature of the nonlocal terms.
Our main effort in this paper is to demonstrate a detailed development of the gradient damage model developed by Aslan and Forest for single crystals under large deformations [18], in parallel with the thermodynamical framework presented in [17]. For that purpose, first of all, the micromorphic damage model which introduces several crystallographic damage mechanisms is demonstrated in Section 2 and the emphasis is placed on the presentation of thermodynamical consistency of the model and identification of the dissipative or energetic character of the gradient variables which is presented in Section 3. The constitutive equations are presented in detail in Section 4 and a specialized modeling case is scrutinized in Section 5.
The second goal of this paper is to use our theory to report on several numerical simulation outcomes generated from model problems. The proposed theory is implemented for the commercial finite element software Abaqus/Standard by implementing a user element subroutine (UEL). In addition to the standart displacement degrees of freedom, a microdamage variable x d is introduced as an additional degree of freedom. The implementation is based on a standard four-node isoparametric element on the interpolated fields with C 0 continuity. Representative numerical examples that illustrate the capacity of the theory and its implementation to damage softening response which results in intense localization are provided together with discussions in Section 6.

Basic Kinematics
Considering a body in the reference configuration represented by X is mapped to the spatial point x at time t, the motion of the body is represented by the smooth function x = χ(X, t). Accordingly, the deformation gradient, the velocity and the velocity gradient are given by Following the The Kröner's decomposition of the deformation gradient [20], the plastic deformation and damage in a crystal is represented in the inelastic distortion and accordingly the decomposition takes the following form: where F e represents the rotation and stretch as elastic distortion and F i is the inelastic distortion capturing the plasticity and damage in the crystal. Figure 1 represents such a deformation for initial and final configurations. However, for the sake of simplicity and the clarity of the proposed theory, the plastic part of the inelastic distortion is neglected throughout this work. Then the deformation gradient simply reads F d represents the local inelastic deformation due to damage generated in the material and invariant under a change in frame. In this framework, the velocity gradient can be decomposed in the following form: where volumetric contribution to elastic and inelastic deformation is defined as Considering the right polar decomposition of the elastic distortion in terms of the right stretch tensor U e and the rotation tensor R e . The right elastic Cauchy-Green tensor and the stretch tensor are defined as Hencky strain [21] is chosen as a strain measure due to its good agreement with the mechanical response of a wide class of materials [22,23].
The elastic and inelastic spin and stretching tensors are Assuming irrotational damage evolution, i.e., The inelastic deformation gradient rate becomeṡ The damage stretching tensor is defined as whereδ s is the damage rates initiated at each damage system which is a scalar internal variable and N d s is the damage flow direction tensor. Finally, the scalar damage rate takes the following forṁ Inspiring from the work of Forest (2009), χ d is defined as the micromorphic counterpart for the damage variable d, for the purpose of mathematical regularization as an additional kinematical degree of freedom. Note that χ d is a positive valued scalar variable which constitutes a subset of micromorphic continuum specifically termed as microdamage continuum by Aslan and Forest (2009). Moreover, the calculation of its gradient ∇ χ d numerically comes at ease, since χ d is defined as a degree of freedom. In order to further develop our theory based on the principle of virtual power, the power expended over the rate of each microvariable is defined as follows: The power overḋ is expended by the microscopic stress π, χḋ is expended by the scalar microscopic stress p, ∇ χḋ is expended by the microscopic stress vector b, and χḋ is expended by the microscopic traction a at the boundary.

Principle of Virtual Power
In light of the definitions provided in Section 2.2, neglecting the inertial terms and assuming that the variation of any independent variable will cause an energy change in the system, external and internal power densities for any given part P are defined as follows where and T R is the 1st Piola-Kirchhoff stress tensor. Principle of virtual power dictates the balance of external and internal powers for any given part P and for all generalized virtual rates consideringḋ ≡ 0 and χḋ ≡ 0 and applying divergence theorem leads to the macroscopic force balance and the traction condition Div Consideringχ ≡ 0 andḋ ≡ 0, choosing an arbitrary χḋ field and applying divergence theorem leads to the microscopic force balance and the traction condition Following the work of Anand et al. [19], local imbalance of energy is formalized aṡ where ψ stands for the free energy per unit volume and C e is introduced into the Quation (24) aṡ where Finally, dissipation density reads

Constitutive Equations
Before the presentation of the constitutive equations, several stress measures to be used in the model need to be introduced.
The Mandel stress is defined by Accordingly, the second Piola stress reads Then, the Cauchy stress takes the following form The first Piola stress is expressed as The free energy is primarily defined as: and from the normality, the energetic state relations can be derived as Introducing the state relations into the dissipation density D defined by Equation (27), the dissipation inequality reads: Therefore, the dissipative terms can simply be defined as: Finally, the dissipation inequality takes its final form as:

Specialization of the Constitutive Equations
Considering the crystallographic structure of a single crystal, the damage evolution is formalized as an addition of each damage system such as: For three mode of fracture, damage evolution is calculated separately asδ s cδ s 1 ,δ s 2 and the number of damage planes N d planes is determined from a given crystal structure. The opening δ s of crystallographic planes with the normal vector n s can be defined as cleavage damage in the system and once cleavage has started, other damage systems are allowed to initiate for the inplane accommodation along orthogonal directions l s 1 and l s 2 , (Figure 2). Three damage criteria are associated to one opening and two accommodation systems: where Y is the damage threshold and the scalar damage is simply the accumulation of damage generated in all systems.ḋ = |δ s c | + |δ s 1 | + |δ s 2 | The correspond equivalent stresses projected on to the damage system takes the form: N s c = n s · M e · n s , N s i = n s · M e · l s (45) For this specific model, the free energy function is taken as a quadratic potential of elastic strain, E e , damage, d, microdamage, χ d and its gradient, ∇ χ d where µ, κ, B and β are positive material constants. Assuming p dis = b dis = 0 Considering the damage criterion: microscopic stress reads the following: Recalling the constitutive relation for the microstress π = π dis + π en , the damage threshold function also becomes: Further, recalling Equation (22), the microscopic force balance takes the from of Helmholtz equation which is postulated as implicit gradient theory for plasticity and damage [11,13,14], where the generalized stresses p and b are not explicitly introduced and the microvariables are defined as nonlocal variables. Div which yields into the Helmholtz equation: Note that the term 1 2 B(d − χ d) 2 in Equation (46) introduces an energetic coupling between d and its micromorphic counterpart χ d, as it is performed in the classical micromorphic theory. Such a penalization necessitates a relatively high coupling modulus, B and the same energetic term also appears in Equation (48) which also serves as a flow rule for d. In the classical micromorphic theory such an energetic coupling is necessary since the elastic strain tensor is coupled with its higher order counterpart and the enrichment is always active even though the behavior is elastic. However, if the same procedure is applied to an internal variable which is dissipative in nature like d, an energetic penalization becomes unnecessary and also brings certain complexity, especially for softening materials. A bulky energetic term in the microforce balance like B(d − χ d) 2 ) cannot contribute to the total dissipation; therefore, plastic softening is always limited by the stored energy through that term. Moreover, for an anisothermal theory, such an energy stored in the body cannot be clearly associated to a physical process in a standard material. In order to tackle that problem, an alternative set of moduli need to be introduced.
In the alternative approach the penalization is not energetic; therefore, coupling modulus, B, does not take place in the free energy function. However, the quadratic coupling between the scalar microvariable and its macroscopic counterpart is preserved: where I is a unit elastic modulus which converts unit strain into unit energy and β * is equivalent to β/B of the conventional approach. Note that in this framework the higher order terms 1 2 I(d − χ d) 2 and 1 2 β * |∇ χ d| 2 are very small compared to the standard elasticity terms, thus, the energy stored through higher order continuum can be ignored. Accordingly, the energetic terms become: and the generalized balance law takes the form: An important observation from Equation (57) is the Laplacian form constructed in the equation provides regularized solutions only for the microvariable χ d and d acts only as a source term. Therefore, in order to attain a similar solution on d, a dissipative penalization is introduced in the microforce balance, such as:σ Considering the stress units in Pascal, for the majority of mechanical problems, Y(d) is in the order of MPa and the B is a penalization constant in the order of elastic moduli. For a high value of B, a satisfactory solution for the local integration is possible only if d follows χ d, otherwise the iterative values ofσ overshoots the elastic trial stressesσ tr and the microforce balance will be violated. In that sense B is selected according to the desired difference between d and χ d. Therefore, the selection criteria for the parameter B is more clear in that context and the energetic term I(d − χ d) can always be ignored since B I.

Numerical Results and Discussion
The model is implemented into a special user element subroutine so-called UEL in ABAQUS software, where the software is used solely as a solver and a visualizer. For the plane strain analysis, 2D four noded quadrilateral elements with bilinear shape functions are used. Microdamage variable χ d is implemented as an extra degree of freedom, in addition to the standard displacements (u 1 and u 2 ). The parameters used in the simulations are provided in Table 1.

Analysis of a 1D Bar
The regularization of a localized deformation with the proposed approach is demonstrated in a numerical model problem, where a bar with a length of 2.5 mm under tension is studied (Figure 3). An opening plane, normal to the direction of the bar is defined and at the center element, an initial defect is introduced where the damage threshold value is slightly reduced. Considering a rate independent damage evolution with linear isotropic softening and only active for opening mode, a microforce balance and the generalized balance law for the alternative case take the following form: whereσ is equivalent tensile stress, H is the softening modulus and S a is the initial material resistance. Using Equation (52) one can rewrite Equation (59) as follows: Further simplifying Equation (60) where solving Equation (61) for χ d gives: Note that for a 1D problem, Equation (63) is a second order homogeneous differential equation with a general form: where which provides two types of solution depending on the sign of ω Note that for a softening problem (providing I > |H|) the solution is sinusoidal an the wave length reads 2π The numerical result presented in Figure 3 is well in line with the analytical solution provided in Equation (68).

Analysis of a 2D Single Crystal Block
As a 2D example, a plate representing a 2D single crystal block under uniaxial tension with a horizontal cleavage plane is investigated (Figure 4). The localization is triggered with an initial geometric defect on the left edge and the cleavage plane is kept on the horizontal axis. The wave length of the localisation band is manipulated by the varying values β. FEA results show that localization path is well in line with the cleavage plane and the localized deformation is captured with in a band whose size is mainly controlled by the length scale l defined in Equation (68).   In order to demonstrate the regularized damage accumulation on a rotated plane, a 2D block with a cleavage plane tilted 25 degrees from the horizontal axis is analyzed under tension with a central defect. Analysis results show that the localization band is perfectly placed on the predefined rotated cleavage plane and the size of the band is kept constant throughout the crack as expected ( Figure 6).

Concluding Remarks
In this work, we have made an attempt to unify the thermodynamically consistent representation of the micromorphic theory presented in the previous work [17] and the microdamage approach proposed by Aslan and Forest [18] for the purpose of modeling crack growth and damage regularization in crystalline solids. It has been shown that the proposed approach is suitable for large deformations and can conveniently be implemented for finite element analysis.
Our analytical effort also demonstrates that the characteristic solution which captures the intrinsic length scale for the localization band can also be achieved by introducing a negligible energetic term. That modification enables the material to undergo a fully damaged state where the equivalent stress reads zero and naturally eliminates the problem of energetic resistance against softening, observed in conventional strain gradient approaches.