A model for the interaction of dislocations with planar defects based on Allen-Cahn type microstructure evolution coupled to strain gradient elasticity

In classical elasticity theory the stress-field of a dislocation is characterized by a $1/r$-type singularity. When such a dislocation is considered together with an Allen-Cahn-type phase-field description for microstructure evolution this leads to singular driving forces for the order parameter, resulting in non-physical (and discretization-dependent) predictions for the interaction between dislocations and phase-, twin- or grain-boundaries. We introduce a framework based on first strain gradient elasticity to regularize the dislocation core. It is shown that the use of strain energy density that is quadratic in the gradient of elastic deformation results in non-singular stresses but may result in singular driving forces, whereas a strain energy, which is quadratic in the gradient of the full deformation tensor, regularizes both stresses and driving forces for the order parameter and is therefore a suitable choice. The applicability of the framework is demonstrated using a comprehensive example.


Introduction
Phase field approaches have proven to be very powerful for the investigation of the formation and evolution of microstructures due to solid-solid phase transformations and twinning. This appears to be the natural framework for the investigation of the interaction of planar crystal defects such as phase-or twin-boundaries with line defects (dislocations, disclinations). A typical phase field model for diffusionless (martensitic) transformations comprises of evolution equations of Allen-Cahn-type for the order parameters φ β where M and α are constants, ρ denotes the mass density, and ψ is a bulk specific free energy. The subscript β indicates the number of the phase, grain or twin variant. Assuming a small perturbation setting, the linear strain tensor 1 E can be additively decomposed into elastic E e and inelastic (i.e., eigenstrain) E in (φ β ) contributions, such that E e E , φ β = E − E in (φ β ). Assuming linear elasticity, the stress S is given by S = C : E e E , φ β , and the specific free energy takes the form ψ E , φ β , θ = 1 2 E e E , φ β : C : E e E , φ β + ψ b φ β , θ .
As a consequence, the evolution equation (1) can be rewritten as In linear elastic Volterra theory, the stresses diverge as the dislocation line is approached. In particular for dislocations the singularity is of 1/r-type. As per Eq.
(3), this results in singular driving forces for the evolution of the order parameters, effectively negating the concepts such as a nucleation barrier or a pile-up stress. Different approaches to regularize the stress in the core region exist in literature based either on the concept of a distributed Burger's vector (Lothe, 1992;Cai et al., 2006), which are inspired by richer microscopic models for dislocations (Peierls, 1940;Nabarro, 1947), or generalized continuum theories (Lazar et al., 2005(Lazar et al., , 2006Lazar and Po, 2015;Po et al., 2018). However, the first strain gradient approach advocated by Po et al. (2018) has the advantage that the obtained regularization is independent of the type of defect in question and therefore does not require any defect-specific information for the determination of model parameters. In principle, these parameters can directly be obtained from atomistic interaction potentials (Admal et al., 2017).
The purpose of this work is to follow a micromorphic approach and to derive a framework which consistently couples first strain gradient elasticity to Allen-Cahn-type microstructure evolution ensuring non-singular driving forces on the order parameters in the presence of line defects.

Balance equations and boundary conditions
The principle of virtual power (PVP) provides a systematic way of deriving field equations and boundary conditions for arbitrary mechanical and coupled problems (cf. Maugin, 1980;Germain, 1973;Del Piero, 2009). In the present work it is used in the following form: The virtual power of the inertia forces P * a balances the virtual power P * int of the internal and P * ext of the external forces acting on any sub-domain S of the material body B for any admissible virtual velocity field v * and virtual rate of order parameter fieldφ * β , i.e., A : B := trA · B , where B is the transpose of B and tr(·) denotes the trace operator. Similarly, we denote third order tensors A by bold calligraphic capital letters and " . . . " is the corresponding scalar product. We use black-board capital letters C for fourth-order tensors. Whenever index-notation is used, summation over latin indices appearing twice is implied and spatial derivatives are denoted using the comma operator, e.g.
For the sake of simplicity we disregard any higher order inertia terms Mindlin (1964) as well as inertial forces acting on the order parameter, resulting in The power of internal forces is given by with L * := grad v * . Here S and T are the Cauchy and higher order stresses, respectively, while π β and ξ β are thermodynamic forces that directly correspond to the internal microforce and microstress introduced by Gurtin (1996). We note that the invariance requirement of P * int with respect to superimposed rigid body motions is satisfied sufficiently by assuming S = S and T · a = (T · a) for arbitrary vectors a. For the power of external forces we consider the very simple case of no body or contact forces acting on L * and gradφ * β , and only a contact (micro)force ζ β actingφ * In order to obtain the consequences of the PVP, the integrals in Eq. (6) are transformed using the following identities and the divergence theorem, resulting in Introducing the surface gradient operator where ∂ n is the directional derivative in the direction of the outward normal n, the third integral in expression (12) can be rewritten as Finally, applying the surface divergence theorem and, for the sake of simplicity, neglecting any wedge line and corner contributions, we find Enforcing Eq. (4) we arrive after a number of straight forward algebraic manipulations at the following field equations and boundary conditions on ∂B We note that, introducing the total stress the balance of linear momentum (16a) regains its standard form for simple materials which is convenient for the numerical implementation.

Constitutive equations
The following equations are formulated assuming a geometrically linear setting, i.e., the displacement gradient is considered to be small || grad u|| 1. In this case the deformation is characterized by the linear strain tensor E = 1 2 grad u + (grad u) . Its gradient will be denoted by Y := grad E.

Laws of state
We choose the following ansatz for the specific free energy and thermodynamic forces The second law of the thermodynamics in the form of the Clausius-Duhem inequality given for the isothermal case by can be exploited using the classical Coleman-Noll procedure to arrive at the laws of state and the residual dissipation inequality 3.2. Free energy and dissipation potential As customary in phase field models for solid-solid transformations, the specific free energy can be split into an elastic, a bulk chemical and an interface contribution as indicated by the subscripts "e" (elastic), "b" (bulk chemical) and "i" (interface). In our formulation, the elastic free energy is of Helmholtz-type, i.e., is the elastic strain, C(φ β ) the stiffness tensor and Λ(φ β ) a gradient length scale tensor (cf. Po et al., 2018). The specific choice of functional dependence of E in (φ β ), ψ b φ β , θ and ψ i φ β , grad φ β , θ on the order parameter φ β is of no relevance at this point; however, we will assume that the interface energy is of the form Using the laws of state (20) we immediately find and combining the first two equations Equation (17) can now be used in two ways: In conjunction with the laws of state (26a) and (26b) it is a constitutive equation for the total stress S t , which enters the balance of linear momentum (18) When combined with Eq. (27), Eq. (17) can be used to determine the true stress S from the total stress S t In order to complete the phase field formulation we require a constitutive equation for π d β , which is obtained in the spirit of classical irreversible thermodynamics aṡ from a dissipation potential Ω π d β that is homogeneous of degree two where M is the so called mobility constant. Combining equations (16b), (21), (26c), (33) and (34) we find the classical Allen-Cahn equation or, explicitely writing down the partial derivatives of ψ, Note that all terms that appear in the driving force, and as per Lazar et al. (2005) the Cauchy stress S in particular, are non-singular even in the presence of dislocations. Interestingly, this is not true for an elastic specific free energy that

Formulation for specific cases
For phase transformations the crystal lattice on both sides of the interface will, in general, be different leading to different elastic properties and a different shape of the dislocation core. In this case the equations (18), (29), (31) and (36) retain their full complexity. However, the strength of the general formulation is that it also covers simplified special cases. In the following, we consider scenarios for which these equations can be strongly reduced and which therefore elucidates the structure of the whole formalism.

Homogeneous bulk material
In the bulk phase the order parameter does not vary in space, i.e., grad The Allen-Cahn equation is fulfilled automatically and Eqs. (31) and (29) For materials with cubic symmetry the gradient length scale tensor Λ is isotropic, i.e., Λ = l 2 I, and the above expressions can be further simplified to the form derived by Lazar et al. (2005)
In the absence of inelastic strain, we have E in (φ β ) = 0. For this case Eqs. (31), (29) and (36) take the form with and The isotropy of the gradient length scale tensor Λ for cubic crystals implies that Q(φ β ) * Λ = Λ = l 2 I, which simplifies Eqs. (38) to the following form and

Twin boundaries and boundaries between grains with inelastic strain
Since the twin variants on both sides of the boundary are related by mirror and/or rotational symmetry transformations between the unit cells, we can -as in the case of grain boundaries -assume that the bulk chemical energy remains unchanged, i.e., ψ b φ β , θ = ψ b (θ), and the elastic stiffness C(φ β ) and the gradient length scale tensor Λ(φ β ) can be expressed using an orthogonal tensor Q(φ β ) as C(φ β ) = Q(φ β ) * C and Λ(φ β ) = Q(φ β ) * Λ, respectively. Under these assumptions we find and For cubic lattices these expressions simplify to and

Phase boundaries between cubic phases
In the case of phase boundaries between different cubic phases the gradient length scale tensor Λ is isotropic on both sides of the interface, even though not necessarily constant across the interface, i.e., Λ = l(φ β ) 2 I. This allows us to reduce Eqs. (31), (29) and (36) to the following form and

Examples
To demonstrate the key properties of the above model numerical simulations using the finite element method are performed using the commercial software "COMSOL Multiphysics" 2 . A uniform mesh with quadratic 3 , quadrilateral elements is used for the domain discretization. The element size is 0.2 nm. Time stepping is performed using the BDF method. Based on the assumptions of the small perturbation hypothesis 4 (Maugin, 1992), we apply traction boundary conditions to the undeformed geometry whenever required. We assume elastostatics with an isotropic stiffness tensor

Regularization in the dislocation core
As shown in Sec. 3.3.1, the present model reduces to the set of equations proposed by Po et al. (2018) in the homogeneous bulk phase. Here, we apply this formulation to a single edge dislocation in an infinite elastic medium: Fig. 1 shows the shear stress component S 12 in the plane perpendicular of this dislocation with and without regularization (l = 2 Å). In the "classical" case without regularization, the stress in the dislocation core is singular, whereas it is well defined and finite for the regularized solution, in analogy to what one would expect from a real atomistic configuration.
2 https://www.comsol.com/ 3 Independent of the chosen shape functions, Comsol does not provide third spatial derivatives of the degrees of freedom. Therefore, in order to obtain the second derivative of strain (third spatial derivative of the displacement), the "Distributed ODE" feature is used in order to introduce additional degrees of freedom, corresponding to the second spatial derivatives of the displacement. For this "Distributed ODE" linear shape functions are employed. 4 Both the displacement u as well as the displacement gradient are considered to be small, i.e., |u| L and || grad u|| 1.

Effect of the regularization on the interaction of dislocations with a moving interface
This examples demonstrates the interaction of dislocations with a moving interface between phase variant 1 (indicated by the superscript "V1") and variant 2 (indicated by the superscript "V2"). The phase mesostructure is described by one single order parameter φ. The only difference between the two variants is with respect to the eigenstrain induced by the phase transformation. This inelastic strain is given as a function of the order parameter by where E in (φ β ) V1 and E in (φ β ) V2 are the eigenstrains of the phases V1 and V2, respectively, ϕ(φ) is a polynomial chosen in accordance with Levitas and Preston (2002) The symmetric bulk chemical free energy takes the following form and the interface energy density is assumed as  For this specific case the resulting set of partial differential equations (18) and (42a) can be further simplified to and These equations are solved for the displacement field u, the order parameter φ and the true stress S. All parameters and coefficients occurring in the above equations are summarized in Tab. 1. The resulting interface energy, computed for a stationary flat interface, is γ = 0.22 J/m 2 . The timescale in the simulation is controlled by the mobility constant M.
Since the simulation time can be arbitrarily re-scaled using the mobility, in our simulations we treat it as dimensionless pseudo time.
The following scenario considers an initially flat interface between variants V1 and V2, and a periodic, immobile dislocation structure with a dislocation spacing of 10 nm within variant 2. For a pictorial representation, see Fig. 2a.
The structure is assumed to be infinite in vertical direction, allowing us to reduce the simulation domain to the dashed 40 nm wide and 10 nm high box in Fig. 2a  peak stresses in the dislocation core (see Fig. 2b). Fig. 3 shows the positions of the V1-V2 interface for different points in simulation time. As the interface approaches the dislocations it bows out due to the interaction with the stress field of the dislocation core. The smaller regularization length results in a larger stress magnitude in the dislocation core region, leading to an arrest of the interface (see Fig. 3a). For the larger regularization length the stress in the vicinity of the dislocation is low enough in order to allow the interface to pass over the dislocation as shown in Fig. 3b.
To analyze the temporal evolution in more detail Fig. 4 visualizes a number of different aspects of the investigated system. The overall V1 phase content, i.e., the area containing phase variant V1 divided by the whole area, as a function of (pseudo) time is shown in Fig. 4a for the two different regularization lengths. There, the most obvious characteristic is the arrest of the interface for both values of l happening simultaneously shortly after t = 2 µs. While the system with the smaller regularization length has already reached a stationary state, the other system shows that the interface "detaches" from the dislocation and swipes the same area per time as before, which shows in the same inclination of the respective line in Fig. 4a.
How does the rate of the V1 phase content evolution change shortly before and after the arresting of the interface took place? In Fig. 4a the peaks at t ≈ 2 and ≈ 3 µs indicate that the phase boundary is accelerating towards the dislocation until its velocity is significantly reduced in the vicinity of the dislocation. The second peak shows that the interface effectively accelerates again after passing the dislocation. This stage is followed by another dip (in between ≈ 3 and ≈ 3.5 µs) where the interface motion right of the dislocation is again decelerated. This behaviors is also  visualized in Fig. 4c, which shows the interface position at equidistant points in time.
Two different phenomenons operate here, which can be understood from the change of sign of the shear stress field of an edge dislocation as shown in Fig. 2a. Recall that the inelastic strain of the interface is governed only by the shear components of the strain tensor. Once the phase boundary is getting close enough to interact with the dislocation, the upper and lower sections of phase V2 are in the regions 3/6 of the dislocation (compare Fig. 2a). The driving force is effectively directed in positive x-direction and causes the acceleration. The central regions of the phase boundary, is located in region 4 of the dislocation stress field and therefore experiences a net driving force that is directed in This is further shown in Figs. 4d and 4e, which relate the motion of the curved interface to the motion of a flat interface for the case without a dislocation structure (indicated by the dotted line).

Summary
In this paper we developed a framework for coupling a phase-field description of planar defects such as phase or twin-boundaries with a discrete representation of dislocations within (anisotropic) first-strain-gradient elasticity. Its main features and advantages in contrast to phase-field within classical elasticity are: • Non-singular stresses at the dislocation core that can be easily calibrated to match molecular statics predictions using the approach of Admal et al. (2017) • Non-singular driving forces for the evolution of the phase-field evolution in the presence of dislocation. This ensures a mesh-independent numerical solution and is a necessary condition for modeling the interaction of dislocations with interfaces such as phase-, grain-or twin-boundaries.
We have shown that in order to ensure regularized driving forces in the dislocation core, a Helmholtz-type elastic free energy that is quadratic in the gradient of the total rather than the elastic strain must be used.
We implemented the proposed framework in the Comsol Multiphics Modeling Software and demonstrated its feasibility and basic properties based on a number of examples. Coupled to a dislocation-dynamics code, we expect this phase-field framework to be a valuable tool for understanding microstructure-evolution on a small scale. Finally, we find the expression where the total stress S t appears in the driving force. In general, this stress cannot be assumed to be bounded in the dislocation-core. This is illustrated in Fig. A.5 that shows the maximum shear stress in the dislocation core for different "thicknesses" of the dislocation, i.e., different discretizations. While the true stress S 12 does not change noticeably once the discretization is sufficiently fine, the total stress S t12 keeps increasing with decreasing thickness of the dislocation.