Mesoscale model of dislocation motion and crystal plasticity
Audun Skaugen Luiza Angheluta Jorge Viñals




Njord Center, Department of Physics, University of Oslo, P.O. 1048 Blindern, 0316 Oslo, Norway

School of Physics and Astronomy, University of Minnesota, 116 Church St. SE, Minneapolis, MN 55455, USA


A consistent, small scale description of plastic motion in a crystalline solid is presented based on a phase field description. By allowing for independent mass motion given by the phase field, and lattice distortion, the solid can remain in mechanical equilibrium on the timescale of plastic motion. Singular (incompatible) strains are determined by the phase field, to which smooth distortions are added to satisfy mechanical equilibrium. A numerical implementation of the model is presented, and used to study a benchmark problem: the motion of an edge dislocation dipole in a hexagonal lattice. The time dependence of the dipole separation agrees with classical elasticity without any adjustable parameters.



46.05.+b,61.72.Bb,61.72.Lk,62.20.F-

Mesoscale model of dislocation motion and crystal plasticity
    Jorge Viñals
    December 30, 2023
============================================================



A phase field theory of a crystalline solid is one of the contending approaches to model defect mediated plastic motion at the nanoscale <cit.>. Due to the diffusive nature of the evolution of the phase field, existing formulations do not maintain mechanical equilibrium during defect motion, and hence are not adequate models of plastic motion. We reinterpret the phase field as the source for singular strains in defected configurations, and add smooth distortions to maintain mechanical equilibrium at all times. The relative motion of two edge dislocations of opposite Burger's vectors under each other's stress field is studied in a hexagonal lattice. In contrast with a direct solution of the phase field equations, we recover the classical law in which the dipole separation scales with time as √(t_0 - t).

Crystal plasticity theory assumes coarse grained volumes that contain a large number of defects, and hence is valid at scales on the order of microns or larger. On the other hand, recent theoretical efforts <cit.> focus on the nanoscale, as new high resolution experiments and large scale simulations are producing a wealth of information about defect motion at the atomic scale. State of the art Bragg Coherent Diffractive Imaging can determine atomic scale displacements with ≤ 30 nm resolution <cit.>. Advanced image processing methods allow the determination of the strain field phase around a single defect, clearly evidencing its multivalued nature. Indeed, single dislocations have been successfully imaged and their motion tracked quantitatively just recently <cit.>. At the same time, nanoscale experiments are revealing new phenomenology that shows that plastic distortion is a rather complex process, characterized by strain bursts and dislocation avalanches <cit.>. These phenomena lie squarely outside of classical crystal plasticity theory. 


A continuum theory of plasticity starts from the statement of incompatibility of the deformation gradient tensor

    ϵ_ilm∂_l w_mk = α_ik,   
    w_mk = ∂_m u_k

where ϵ_ilm is the anti symmetric Levi-Civita tensor, α_ik the dislocation density tensor, and w_mk the distortion tensor <cit.>. The integral of α_ik over a surface is the sum of the Burger's vectors 𝐛 corresponding to all the dislocation lines n that pierce the surface ∫_Sα_ijdS_j = ∑_n b_i^n. For any given distribution of topological defects in the material, α is fixed, but not the distortion, which can be decomposed into a singular part, the curl of which yields α_ij, and a smooth strain which we denote by u_ij^δ. The smooth strain is compatible ϵ_ikmϵ_jln∂_klu_mn^δ = 0. Regardless of the state of distortion, plastic motion is slow on the scale of lattice vibration, and occurs in mechanical equilibrium, as the stress σ_ij adiabatically follows the instantaneous distribution of dislocations, ∂_jσ_ij = 0. Closure generally requires a constitutive relation involving the stress and the smooth deformation. These considerations and appropriate boundary conditions are sufficient to specify the static problem. Dynamically, over the time scale appropriate for plastic flow, an evolution equation needs to be introduced for the dislocation density tensor <cit.>. In field dislocation dynamics theories, its evolution is kinematically related to the velocity of the dislocation lines, which in turn require a constitutive definition in terms of a local free energy and a dissipation function <cit.>. The phase field crystal model of defect motion, as currently formulated <cit.>, can be used to specify most of the static and dynamic features just described, but not all, as discussed below.

The phase field, a scalar function of space and time, is a physical order parameter that describes the dimensionless mass density of the crystalline phase ψ(𝐫) = ψ_0 + ∑_𝐠 A_𝐠 e^i 𝐠·𝐫, where the sum extends over all reciprocal lattice vectors 𝐠 of the lattice. A non convex free energy functional for an isothermal system ℱ is introduced so that its minimizer ψ^* has the desired symmetry of the crystalline phase. Lattice constants appear as parameters. In dimensionless units, we use ℱ[ψ] = ∫ d 𝐫 f(ψ,∇^2ψ) with  f(ψ,∇^2ψ)= (ℒψ)^2/2+ r^2ψ^2/2+ψ^4/4, and ℒ = 1 + ∇^2 <cit.>. The only remaining constant parameter r is the dimensionless distance away from the symmetry breaking bifurcation. For r > 0, ψ^* = 0 is the only stable solution. For r < 0, and depending on the conserved spatial average ψ_0, ψ^* is periodic with wavenumber unity in our dimensionless units, but of various symmetries. For simplicity, we consider here a two dimensional system where the equilibrium configuration is a hexagonal phase with lattice constant a = 4π/√(3). The Burger's vector density in 2D is B_k (𝐫) = α_3,k(𝐫), k=1,2. Our results, however, can be readily extended to three dimensions. The temporal evolution of the phase field ψ is defined to be relaxational and driven by free energy reduction 

    ∂_tψ(𝐫,t) = ∇^2δℱ/δψ(𝐫,t)

where δ / δψ(𝐫,t) stands for the variational derivative with respect to the phase field. 

For smooth distortions of ψ^*, the free energy ℱ suffices to determine the stress-strain relation <cit.>. For small distortions, we define a non-singular stress σ^ψ <cit.>


    σ_ij^ψ = σ̃^ψ_ij_c,   σ̃_ij^ψ = [∂_iℒψ]∂_jψ - [ℒψ]∂_ijψ + f δ_ij,


with σ̃^ψ_ij given by the local variation of ℱ with ∂_i u_j, and ·_c denoting a spatial average across a region roughly corresponding to a unit cell. σ_ij^ψ is symmetric and related to the strain field u_ij= (∂_i u_j+∂_j u_i)/2 according to linear elasticity. For the hexagonal phase under discussion, the relation is that of isotropic elasticity


    σ_ij^ψ = λδ_iju_kk + 2μ u_ij


with Lamé coefficients λ = μ = 3 A_0^2 <cit.>. The quantity A_0 is the amplitude of the uniform mode in a multiple scale amplitude expansion of ψ^*.

The phase field model has been used extensively to describe static and dynamic defected configurations. Following early work on dislocation motion and grain boundaries in roll patterns <cit.>, the phase field crystal theory has been used to study dislocation <cit.> and grain boundary motion <cit.>. Strain fields have been explicitly extracted <cit.>, or imposed to analyze strained film epitaxy <cit.>, and considered as the limiting case of phonon degrees of freedom <cit.>. More complex properties of defect motion such as specification of slip systems, defect mobilities, and Peierls barriers are also given by phase field kinetics <cit.> thus opening the door to the study of defect pinning, bursts, and avalanches.

However, whereas for a specified and fixed defected configuration minimizers of ℱ with appropriate boundary conditions can be found that are in mechanical equilibrium, any local deformation of ψ(𝐫,t) propagates only diffusively according to Eq. (<ref>). The relevant transverse diffusion constant is small, and can even vanish <cit.>. This is not physical for a crystalline solid, as has been already recognized <cit.>. In ordinary crystals, unlike the phase field model, elastic equilibrium compatible with a transient distribution α_ik(𝐫,t) and boundary conditions is established quickly, in a time scale determined by damping of elastic waves in the medium.

To overcome this difficulty, we propose to use the phase field ψ(𝐫,t) only as an indicator function of defect location and topology, as well as governing local relaxation near defect cores. The field ψ(𝐫,t) determines the source for lattice incompatibility in Eq. (<ref>), the solution of which is only a particular singular solution for the deformation field. A smooth distortion 𝐮^δ (in the null space of the curl) must be added to this particular solution to enforce elastic equilibrium. Equation (<ref>) for the newly displaced phase field ψ^'(𝐫) = ψ(𝐫 + 𝐮^δ) provides for defect motion in a manner that is consistent with the Peach-Kohler force <cit.>. Plastic motion is uniquely specified, with the only constitutive input being the free energy functional ℱ. We discuss in what follows the details of our computational implementation, and specifically address the relative motion of a dislocation dipole in a 2D hexagonal phase. 

We decompose the stress field into a singular part arising from the current phase field configuration, σ^ψ, and a small contribution arising from the smooth distortion σ^δ, so that σ = σ^ψ+σ^δ is in mechanical equilibrium ∇·σ = 0. This condition is satisfied by introducing the Airy function χ, which in two dimensions reads σ_ij  = ϵ_ikϵ_jl∂_klχ. Inverting Eq. (<ref>), we have in 2D

    u_ij = 1/2 (∂_i u_j + ∂_j u_i)
    = 1/2μ(σ_ij - κδ_ijσ_kk),

where κ = λ/2(λ+μ). Inserting Eq. (<ref>) into the incompatibility relation in 2D ϵ_ikϵ_jl∂_klu_ij=ϵ_ij∂_i B_j(𝐫) (e.g. <cit.>) and expressing the stress in terms of χ-function, we obtain that 


    1-κ/2μ∇^4χ = ϵ_ij∂_i B_j(𝐫),

where B_j(𝐫) = ∑_α b^n_jδ(𝐫 - 𝐫_n) is the dislocation density in 2D for a configuration of dislocations with Burgers vector 𝐛^n at locations 𝐫_n. In Ref. <cit.>, we explicitly computed 𝐁(𝐫) through complex demodulation of the phase field ψ(𝐫,t). Demodulation yields both the amplitude and phase of the deformation field; the former going to zero at the defect core, the latter undergoing a discontinuity across a line that terminates at the core. Figure <ref>(a) shows a dislocation dipole in a 2D hexagonal lattice, and Figure <ref>(d) the right hand side of Eq. (<ref>)) obtained by demodulation. We proceed differently here and introduce a more efficient numerical procedure that does not require demodulation. The smooth strain u_ij^δ is compatible (ϵ_ikϵ_jl∂_iju_kl^δ = 0) and therefore, with Eq. (<ref>), the corresponding stress satisfies,

    ϵ_ikϵ_jl∂_ij(σ_kl^δ - κδ_klσ_ll^δ) = 0.

We now proceed as if the linear decomposition σ = σ^ψ + σ^δ holds everywhere, including near dislocation cores as defined by the phase field. This is not strictly correct, but this decomposition results in a distortion u_ij^δ which is everywhere compatible. However, as discussed below, the computed stress field σ will be divergence free only away from any defect core. Given this decomposition σ_ij^δ = ϵ_ikϵ_jl∂_ijχ - σ_ij^ψ, we find an analogous result to  Eq. (<ref>),


    (1-κ)∇^4χ = (ϵ_ikϵ_jl∂_ijσ_kl^ψ - κ∇^2σ_kk^ψ).


Note that the stress σ^ψ (as defined in Eq. (<ref>)) is smooth and bounded, so the right-hand side of Eq. (<ref>) can only give a nonsingular approximation to the singular right hand side of Eq. (<ref>). Figure <ref>(c) shows the right-hand side of Eq. (<ref>) obtained numerically for the dislocation dipole which is in good agreement Eq. (<ref>) obtained through demodulation (Fig. <ref>(d)). Both methods act as regularizations of the singular density at defect cores.

From a given phase field configuration at time t, ψ(𝐫,t), we compute σ^ψ from Eq. (<ref>), and then solve Eq. (<ref>) to obtain χ and therefore σ. The difference σ_ij^δ = ϵ_ikϵ_jl∂_ijχ - σ_ij^δ leads to the smooth strain u_ij^δ = 1/2 μ ( σ_ij^δ - νσ_kk^δδ_ij ) which is, by construction, compatible. It can, therefore, be integrated to obtain a compatible deformation 𝐮^δ. The final step in the computation is to redefine the phase field ψ^' (𝐫,t) = ψ(𝐫 + 𝐮^δ,t). 

Both the stress-strain relation and stress superposition only hold far from defect cores. We define the stress of this newly deformed configuration everywhere as

    σ_ij = σ^ψ_ij + σ^δ_ij = σ_ij^ψ + λδ_iju^δ_kk + 2μ u^δ_ij,

which satisfies ∂_jσ_ij = 0 only far from defect cores, not at short distances. This is not a problem as standard diffusive evolution of the phase field suffices to equilibrate the stress near cores in time. We discuss this further below, and in Fig. <ref>.




The integration of the compatible strain u_ij^δ to obtain 𝐮^δ is carried out through a Helmholtz decomposition into curl-free and divergence-free parts u_i^δ = ∂_i V + ϵ_ij∂_j A. Applying the divergence to this expression, one obtains a Poisson equation for the potential V, ∂_i u_i^δ = u_ii^δ = ∇^2 V, which is easily solved by spectral methods. On the other hand, taking the curl we find 
ϵ_ij∂_i u_j^δ = ϵ_ijϵ_jk∂_ikA = -∇^2 A, which is a Poisson equation for A. Unfortunately the source term depends on the antisymmetric part of the smooth deformation gradient, which we do not obtain directly from the elastic stress, as this only depends on the symmetric part. We therefore apply another Laplacian operator to the equation, and use the compatibility relation ϵ_ij∂_iju^δ_k = 0 to find


    ∇^4 A = -ϵ_ij∂_ik(∂_k u_j^δ + ∂_j u_k^δ)
    = -2ϵ_ij∂_iku_jk^δ.


This is a biharmonic equation for A with a known source term, which is again easily solved by spectral methods. In particular, if k_i are the components of the 𝐤 vector and û_ij^δ are the Fourier components of the residual strain, the Fourier components of the residual deformation can be expressed as 

    û_i^δ = - ik_i/k^2û_jj^δ + 2iϵ_ijϵ_rsk_jk_rk_l/k^4û_sl^δ,

with the k=0 component chosen to be zero to avoid rigid body displacements. We then compute the distorted phase field ψ'(𝐫) = ψ(𝐫 + 𝐮^δ) on the original grid 𝐫 by expanding in Taylor series up to fifth order in 𝐮^δ.





We present next the results of a numerical study of a benchmark configuration: the relative motion of two edge dislocations along their glide plane, Fig. <ref>(a). The computational domain is assumed periodic containing 120x120 unit cells with a spatial resolution of a/8 = π/2√(3) in the x direction, and 2π/7 in the y direction. The initial distance between dislocations is 24a, and we consider two sets of values r = -0.2, ψ_0 = 0.265 and r = -0.8,ψ_0 = -0.43. We prepare the initial condition in the same way as in Ref. <cit.>, and numerically solve Eq. (<ref>) using an exponential time differencing method with a time step of Δ t = 0.1 <cit.>. According to linear elasticity theory for an isotropic medium, the shear stress in such a configuration is,

    σ_xy = 2μ(λ+μ)/λ+2μ∑_αb_x^n/2πcosϕ_n cos(2ϕ_n)/|𝐫-𝐫_n|,

where ϕ_n is the azimuth relative to dislocation n. Figures (<ref>)(a,c) show σ^ψ along a line that includes the two dislocation cores at time t = 700 obtained by direct integration of Eq. (<ref>), and compares it to Eq. (<ref>). Divergences in Eq. (<ref>) are regularized by the phase field, and the stress near the cores is relatively well described by σ^ψ. Far from the cores, however, the two stresses show qualitatively different asymptotic dependence. Figures (<ref>)(b,d) show the stress in a configuration in which the smooth distortion (<ref>) has been applied between time steps. The stress is still regularized near defect cores, yet, by construction, agrees with linear elasticity away from them.

Figure <ref> shows the time dependence of the dislocation velocity as a function of dipole separation as given by direct integration of the phase field model, Eq. (<ref>), Fig. <ref>(a), and by our model, Fig. <ref>(b). For reference, we also show the expected result from elasticity theory by using the Peach-Kohler force with stress (<ref>), and mobility derived from ℱ (Eq. (45) in Ref. <cit.>). There are no adjustable parameters in the calculation of the analytic velocity. The two dislocations move towards each other until they annihilate, with a velocity inversely proportional to their separation. Our model captures this result well for a range of parameters r, with slight stick-slip motion observed at larger |r|, visible as oscillations in the dislocation velocity.





To summarize, we have argued that the phase field crystal model currently in use lacks deformation as an independent variable, and as a consequence fails to maintain proper mechanical equilibrium during plastic motion. We retain the model because it provides for lattice and topological defect structures as derived properties from the phenomenological free energy. It also allows regularization of defect cores and singular stresses. Phase field kinetics is also consistent with the classical Peach-Kohler force, with mobility that is again specified by the free energy ℱ. We take the view, however, that the phase field is not adequate to describe the distortion of the lattice away from moving defect cores, and hence supplement it with a smooth distortion field, compatible with the topological content of the phase field, but defined so as to maintain mechanical equilibrium everywhere away from defect cores. When the evolution of ψ'(𝐫,t) is thus constrained to satisfy mechanical equilibrium at all times, we show numerically that our model agrees with the classical law of motion for a dislocation dipole in isotropic, linear elasticity. These results put the phase field crystal model on firmer ground to study more complex defected configurations at the nanoscale.

This research has been supported by a startup grant from the University of Oslo, and the National Science Foundation under contract DMS 1435372.

apsrev4-1


