Landau theory of crystal plasticity

TL;DR

This paper develops a Landau theory-based nonlinear continuum elasticity model for crystal plasticity, capturing dislocation nucleation and defect complexity through symmetry considerations and athermal dynamics.

Contribution

It introduces a novel Landau theory framework with infinite energy wells based on crystal symmetry for modeling plastic flow in crystals.

Findings

Dislocation nucleation modeled in 2D crystals under shear.

Symmetry group GL(3,Z) influences defect configurations.

Higher symmetry phases impact plastic deformation processes.

Abstract

We show that nonlinear continuum elasticity can be effective in modeling plastic flows in crystals if it is viewed as Landau theory with an infinite number of equivalent energy wells whose configuration is dictated by the symmetry group GL(3,Z). Quasi-static loading can be then handled by athermal dynamics, while lattice based discretization can play the role of regularization. As a proof of principle, we study in this Letter dislocation nucleation in a homogeneously sheared 2D crystal and show that the global tensorial invariance of the elastic energy foments the development of complexity in the configuration of collectively nucleating defects. A crucial role in this process is played by the unstable higher symmetry crystallographic phases, traditionally thought to be unrelated to plastic flow in lower symmetry lattices.

Equations

$$f_d(\tilde{\mathbf{C}})=\beta {\psi }_1(\tilde{\mathbf{C}})+{\psi }_2(\tilde{\mathbf{C}})$$

$$\nabla \mathbf{y}=1+\alpha ({\mathbf{e}}_1\otimes {\text{e}}_1^{\perp }).$$

Full text

Keywords

Crystal plasticity is the simplest among yield phenomena in solids Bonn et al. (2017), and yet it has been compared in complexity to fluid turbulence Cottrell (2002); Choi et al. (2012). The intrinsic irregularity of plastic flow in crystals Vinogradov and Estrin (2018) is due to short and long range interaction of crystal defects (dislocations) Zepeda-Ruiz et al. (2017) dragged by the applied loading through a rugged energy landscape Salman and Truskinovsky (2012a); Kubin (2013); Ispánovity et al. (2014). Fundamental understanding of plastic flow in crystals is crucial for improving hardening properties of materials Hughes and Hansen (2018), extending their fatigue life Irastorza-Landa et al. (2016), controlling their forming at sub-micron scales Csikor et al. (2007) and building new materials McDowell (2018a).

Macroscopic crystal plasticity relies on a phenomenological continuum description of plastic deformation in terms of a finite number of order parameters representing amplitudes of pre-designed mechanisms. These mechanisms are coupled elastically and operate according to friction type dynamics Roters et al. (2010); Forest et al. (2019); Gurtin and Anand (2009); Han et al. (2005); Ste (2018). The alternative microscopic approaches, relying instead on molecular dynamics Juan and Kaxiras (1993); Bulatov et al. (1998); Lyu and Li (2019); McDowell (2018b); Moretti et al. (2011); Skaugen et al. (2018); Tarp et al. (2014); Elder et al. (2002), can handle only macroscopically insignificant time and length scales McDowell (2019). An intermediate discrete dislocation dynamics approach focuses on long range interaction of few dislocations, while their short range interaction is still treated phenomenologically Shilkrot et al. (2004); Song et al. (2019); El-Awady et al. (2016). Collective dynamics of many dislocations can be also described by the dislocation density field, however, rigorous coarse-graining in such strongly interacting system still remains a major challenge El-Azab (2000); Xia and El-Azab (2015); Groma (2019); Acharya and Roy (2006); Chen et al. (2013); Sandfeld et al. (2011); LeSar (2014); Valdenaire et al. (2016).

A highly successful computational bridge between microscopic and macroscopic approaches is provided by the quasi-continuum finite element method which uses adaptive meshing and employs ab initio approaches to guide the constitutive response at different mesh scales Tadmor et al. (1996); Tadmor and Miller (2017); Kochmann and Amelang (2016); Li et al. (2002, 2004). Its drawbacks, however, are spurious effects due to matching of FEM representations at different scales and a high computational cost of reconstructing the constitutive response at the smallest scales Luskin and Ortner (2013).

In this Letter we propose a synthetic approach dealing with the macroscopic quantities such as stresses and strains, while accounting correctly for the exact symmetry of the crystal lattice. Our main assumption is that meso-scale material elements are exposed to the periodic energy landscape which resolves lattice-invariant strains, including shears related to slip Bulatov and Argon (1994); Kaxiras and Boyer (1994), see Fig. 1. Our approach follows the original proposal by Ericksen that the energy periodicity in the space of tensors should be made compatible with geometrically nonlinear kinematics of crystal lattices Ericksen (1970, 1973, 1977, 1980), and we also build upon subsequent important developments of the mathematical formalism in Parry (1976); Folkins (1991); Parry (1998); Pitteri and Zanzotto (2002); Conti and Zanzotto (2004).

This general program can be viewed as far reaching generalization of the Frenkel-Kontorova-Peierls-Nabarro model accounting for energy periodicity along a single slip plane Frenkel and Kontorova (1939); Peierls (1940); Nabarro (1947). Scalar models with periodic energies, dealing with multiple slip planes, have been used before to describe dislocation cores Kovalev et al. (1993); Landau (1994); Carpio and Bonilla (2003), to simulate dislocation nucleation Plans et al. (2007); Bonilla et al. (2007); Lomdahl and Srolovitz (1986); Srolovitz and Lomdahl (1986) and to capture intermittency of plastic flows Salman and Truskinovsky (2011, 2012b). Their tensorial versions with linearized kinematics were considered in Minami and Onuki (2007); Onuki (2003); Carpio and Bonilla (2005); Geslin et al. (2014).

In the proposed kinematically nonlinear theory the role of the order parameter is played by the metric tensor (characterizing local deformation), and the bottoms of the energy wells correspond to lattice invariant deformations. The continuum Landau energy density of this type must be invariant under the infinite symmetry group $\mathrm{GL}(3,\mathbb{Z})$ and since the ground state in this case is necessarily hydrostatic Ericksen (1973); Fonseca (1987), the regularization is necessary. Theories incorporating various elements of the tensorial $\mathrm{GL}(3,\mathbb{Z})$ symmetry have already proved useful in the description of reconstructive phase transitions Dmitriev et al. (1988); Horovitz et al. (1989); Sanati et al. (2001); Bhattacharya et al. (2004); Pérez-Reche et al. (2009); Vattré and Denoual (2016) and in this Letter we extend this idea to the modeling of crystal plasticity proper, see also Biscari et al. (2015).

From the perspective of Landau theory with an infinite number of equivalent energy wells, plastically deformed solid can be viewed as a multi-phase mixture of equivalent phases. Given the large magnitude of the ’transformation strain’, such phases are localized at the scale of the regularizing length and the domain boundaries appear macroscopically as linear defects mimicking dislocations. Plastic yield can be then interpreted as an escape from the reference energy well and plastic mechanisms can be linked to low-barrier valleys in the energy landscape. Friction type dissipation emerges as a result of homogenization of an overdamped athermal dynamics in a rugged energy landscape Puglisi and Truskinovsky (2005); Mielke and Truskinovsky (2011). It is important to notice that such approach incorporates both long and short range dislocation interactions; it also correctly describes plastic slip even though the dislocations cores are regularized and blurred on the scale of the unit cell.

In this Letter we use this general approach to study the peculiarities of collective dislocation nucleation in crystals with different symmetries and show that in low symmetry lattices there is a nontrivial coupling between plastic mechanisms due to the presence of degenerate mountain passes representing lattices with higher symmetry. The crucial role in the formation of self induced plastic disorder is then played by the unstable high symmetry phases, traditionally thought to be unrelated to plastic flow. The general conclusion is that, rather paradoxically, the global crystal  symmetry can induce frustration and become the cause of lattice incompatibility developing during plastic deformation.

Consider a continuous deformation $\mathbf{y}=\mathbf{y}(\mathbf{x})$, where $\mathbf{y}$ are actual and $\mathbf{x}$ are reference coordinates. The energy density of an elastic solid can depend on the deformation gradient $\bm{\nabla}\mathbf{y}$ only through the metric tensor $\mathbf{C}=(\bm{\nabla}\mathbf{y})^{T}\bm{\nabla}\mathbf{y}$. To account for all deformations that map a Bravais lattice into itself we must require that $f(\mathbf{C})=f(\bf m^{T}\mathbf{C}\bf m)$ for any $\bf m$ from a discrete group conjugate to $\mathrm{GL}$ $(3,\mathbb{Z})$ and comprised of all invertible matrices with integral entries and determinant $\pm 1$, see Supplementary Material com . In the presence of such symmetry, the space of metric tensors $\mathbf{C}$ partitions into periodicity domains, each one containing an energy well equivalent to the reference one. If we know the structure of the energy in one of such domains, we can use, for instance, the Lagrange reduction Engel (1986); Conti and Zanzotto (2004) to find its value in any other point.

In the special case of 2D lattices, which we focus on in what follows, a section $\det\mathbf{C}=1$ in the 3D space of tensors $\mathbf{C}$ can be used to visualize the implied tensorial periodicity of the energy, see Fig 2. The global picture is made visible if we map this section into a complex half-plane using the function $z=C_{11}^{-1}(C_{12}+\mathrm{i})$ Folkins (1991); Parry (1998). For instance, the point $\bf S_{1}$ in Fig. 2, corresponding on the complex half plane to $z=\mathrm{i}$ describes a square lattice with the basis vectors aligned with the close-packed directions: ${\bf e}_{1}=(1,0)$, ${\bf e}_{2}=(0,1)$. Simple shear $\mathbb{1}+\,\mathbf{e}_{1}\otimes\textbf{e}_{1}^{\perp}$, where $\textbf{e}_{1}^{\perp}$ is a vector orthogonal to $\textbf{e}_{1}$, maps this point into its symmetric counterpart $\mathrm{i}+1$ (point $\bf S_{2}$); Another square lattice, corresponding to point $\bf S_{3}$ in Fig. 2 with $z=1/2(1+\mathrm{i})$, can be obtained from the lattice $\bf S_{1}$ by the shear $\mathbb{1}+\,\mathbf{e}_{2}\otimes\textbf{e}_{2}^{\perp}$. Instead, the point $\bf T_{1}$ in Fig 2, corresponding to $z=\frac{1}{2}+\frac{\sqrt{3}}{2}\mathrm{i}$, describes a triangular lattice (with hexagonal symmetry) whose basis vectors ${\bf e}_{1}=\gamma(1,0)$, ${\bf e}_{2}=\gamma(1/2,\sqrt{3}/2)$ with $\gamma=(4/3)^{1/4}$ are again aligned with the close-packed directions. Its closest equivalent neighbors are $\bf T_{2}$ and $\bf T_{4}$ corresponding to $z=\frac{1}{2}+\frac{\sqrt{3}}{2}\mathrm{i}\pm 1$. They are reachable from $\bf T_{1}$ by the shear deformation $\mathbb{1}\pm\,\mathbf{e}_{1}\otimes\textbf{e}_{1}^{\perp}$.

To demonstrate the possibility of yield-inducing instabilities in a material with such energy, it is sufficient to consider a system under the most constraining, affine displacement control. Given the gradient nature of the order parameter, the (spinodal) instability of a homogeneous state should be linked to the local loss of rank-one-convexity of the energy Ogden (1997); Grabovsky and Truskinovsky (2014). This is equivalent to the loss of positive definiteness of the acoustic tensor $\mathbf{Q}$ with components $Q_{ik}=A_{ijkl}n_{j}n_{l}$, where $A_{ijkl}=\partial^{2}f/(\partial_{j}y_{i}\partial_{l}y_{k})$ is the fourth-order incremental elastic tensor and $\mathbb{n}$ is a unit vector Van Vliet et al. (2003); Miller and Rodney (2008).

For illustrative purposes we now choose a particular energy density $f=f_{v}+f_{d}$ which decouples into a volumetric $f_{v}(\det\mathbf{C})$ and a deviatoric $f_{d}(\mathbf{C}/(\det\mathbf{C})^{1/2})$ parts. Since $\det\mathbf{C}$ is invariant under $\mathrm{GL}(2,\mathbb{Z})$, our symmetry constraints concern only the deviatoric part $f_{d}$. This function needs to be specified only inside a single periodicity domain with the suitable conditions on its boundary ensuring required smoothness Parry (1976). The lowest order polynomial representation of $f_{d}$, which guarantees the continuity of the elastic moduli, was constructed in Conti and Zanzotto (2004); for the general non-polynomial representation see Folkins (1991).

If the reference lattice is either square or triangular symmetry, the minimal potential can be chosen in the form Conti and Zanzotto (2004):

[TABLE]

where $\psi_{1}={I_{1}}^{4}\,I_{2}-41\,{I_{2}}^{3}/99+7\,I_{1}\,I_{2}\,I_{3}/66+{I_{3}}^{2}/1056,$ and $\psi_{2}=4\,{I_{2}}^{3}/11+{I_{1}}^{3}\,I_{3}-8\,I_{1}\,I_{2}\,I_{3}/11+17\,{I_{3}}^{2}/528.$ The hexagonal invariants here have the structure: $I_{1}=\frac{1}{3}(\tilde{C}_{11}+\tilde{C}_{22}-\tilde{C}_{12})$, $I_{2}=\frac{1}{4}(\tilde{C}_{11}-\tilde{C}_{22})^{2}+\frac{1}{12}(\tilde{C}_{11}+\tilde{C}_{22}-4\tilde{C}_{12})^{2}$, and $I_{3}=(\tilde{C}_{11}-\tilde{C}_{22})^{2}(\tilde{C}_{11}+\tilde{C}_{22}-4\tilde{C}_{12})-\frac{1}{9}(\tilde{C}_{11}+\tilde{C}_{22}-4\tilde{C}_{12})^{3}.$ The choice $\beta=-1/4$ enforces the square symmetry on the reference state, while choosing $\beta=4$ we bias the reference state towards hexagonal symmetry; the energy landscapes in those two cases are illustrated in Fig. 3. The volumetric energy density will be chosen in the simplest form $f_{v}(s)=\mu(s-\log(s))$, which excludes configurations with infinite compression; the coefficient $\mu$ plays the role of a bulk modulus.

The resulting ’yield surfaces’ are shown in Fig. 3. To understand the nature of the associated instabilities, consider the case when a simple shear is imposed on the boundary

[TABLE]

Starting, in the case of square lattice, from the homogeneous reference state $\bf S_{1}$, we find that at instability point the condition $\det{\bf Q}=0$ produces two, almost simultaneously destabilized directions $\mathbb{q}=\bm{\nabla}\mathbf{y}[\mathbb{n}]/|\bm{\nabla}\mathbf{y}[\mathbb{n}]|=(\cos\xi,\sin\xi)$: the first one with $\xi\approx-0.11$rad, almost perpendicular to the deformed $\mathbb{e}_{2}$, and the second one with $\xi\approx 1.55$rad, almost perpendicular the the deformed $\mathbb{e}_{1}$. The near-degeneracy of the bifurcation is an indicator that two ’slip planes’ may be activated. In the case of triangular lattice, the instability along a similar loading path originating at $\bf T_{1}$ produces a single unstable direction $\xi\approx-1.25$rad which is incommensurate with the lattice. In this case one can expect only one ’slip plane’ to be activated. Our numerical experiments show that the acoustic-tensor-based analytical instability conditions are in agreement with direct numerical simulations.

Before addressing the post-bifurcational behavior consider a single edge dislocation trapped by the lattice somewhere far from the boundaries. In Fig. 4 we illustrate the corresponding stress distribution which matches the classical continuum far field with $r^{-1}$ asymptotics while also resolving (at a scale of the mesh) the core region. Solutions like this can be helpful in calibrating the model using molecular statics simulations.

The collective nucleation pattern emerging after a stress drop is illustrated in Fig. 5 for both types of lattices. The results are presented on both, the configurational space, Fig. 5(a,c), where each point corresponds to a single element of the mesh and the actual physical space, Figs. 5(b,d), where the color of the dots (nodal points) indicate the level of stress. The configurational points, all located initially at the bottom of the reference energy well, disperse as a result of the massive nucleation event. The ensuing spatial dislocation distribution is quasi-regular with pile-ups at the rigid boundaries. Note the formation of characteristic entanglements with dislocations on two slip planes blocking each other (in the case of square lattice); there is also some disorder due to unavoidable numerical noise.

Note that in the case of square lattice, the system is driven by the loading device from the reference state $\bf S_{1}$ towards the equivalent state $\bf S_{2}$. At the "yielding" threshold, which marks the end of the elastic regime, the homogeneous configuration $\bf S_{1}$ loses stability and the ensuing pattern represents (outside the core regions) a mixture of three ’pure’ states $\bf S_{1}$, $\bf S_{2}$ and $\bf S_{3}$, see Fig.5(a). While the appearance of the state $\bf S_{2}$ is natural, because the corresponding ’plastic mechanism’ is favored by the loading, the main complexity of the resulting dislocation pattern is due to the emergence of the state $\bf S_{3}$. It indicates the activation of the second plastic mechanism, decoupled (in the nonlinear theory) from the first one.

The appearance of the state $\bf S_{3}$ can be understood if we recall that the linear stability analysis predicted two almost simultaneously unstable modes aligned with the slip directions in the deformed state. While one of these directions is indeed aiming towards the energy well $\bf S_{2}$, the other one, which bifurcates first, is directed towards $\bf S_{3}$. Our numerical simulations show that the latter instability mode grows faster which can be interpreted as, somewhat counter-intuitive, early stage dominance of the secondary ’plastic mechanism’. The flow of configurational points passes near the unstable equilibrium state $\bf T_{1}$, corresponding to a triangular lattice, where it splits into three streams directed towards the configurations $\bf S_{1}$, $\bf S_{2}$ and $\bf S_{3}$, see Movie S1 in com .

This behavior becomes more transparent if we consider a smoother energy potential. Note that the symmetry transformations from $\mathrm{GL}(2,\mathbb{Z})$ correspond on the upper complex half-plane to the fractional (Moebius) transformations with integral entries of the type \big{(}m_{22}z+m_{12}\big{)}\big{/}\big{(}m_{21}z+m_{11}\big{)} Parry (1976); Folkins (1991). This observation links the infinitely periodic energy densities for 2D crystalline materials with the classical modular functions Schoeneberg (2012), with the most well known example provided by the Klein invariant $J(z)$ van Wijngaarden (1953), see Supplementary Material for more details com . One can show that for this holomorphic function $J|_{\bf S_{i}}=1$, $J^{\prime}|_{\bf S_{i}}=0$, while $J|_{\bf T_{i}}=J^{\prime}|_{\bf T_{i}}=J^{\prime\prime}|_{\bf T_{i}}=0$. Therefore, the corresponding potentials with the reference square and triangular lattices can be chosen in the form: $f_{d}(z)=|J(z)-1|$ (square lattice) and $f_{d}(z)=|J(z)|^{2/3}$ (triangular lattice); the exponents are chosen to ensure a non-degenerate linear-elastic response close to the bottoms of the energy wells. The energy landscapes and the yield surfaces for such potentials are qualitatively similar to the ones presented in Fig. 3.

Note that the choice $f_{d}(z)=|J(z)-1|$ for a square lattice turns the ’triangular’ critical point $\bf T_{1}$ and all its symmetric counterparts into degenerate ‘monkey saddles’, characterized by the local Taylor expansion of the form $x^{3}-3xy^{2}$, see Fig.6(a). The flow of configurational points directed initially towards such unstable state (say, $\bf T_{1}$) will therefore necessarily split into three streams directed towards the stable states (say, $\bf S_{1}$, $\bf S_{2}$ and $\bf S_{3}$).

Superficially, the situation looks a bit different in the case of the polynomial energy (1), where the Hessian is nondegenerate at the point $\bf T_{1}$ which corresponds in this case to a shallow energy maximum. However this maximum is surrounded by the three nondegenerate saddles $\bf R_{1}$, $\bf R_{2}$ and $\bf R_{3}$ describing rhombic lattices, see Fig.6(b), and the general conclusion about the activation of the secondary plastic mechanism and the ultimate dispersion over three energy wells $\bf S_{1}$, $\bf S_{2}$ and $\bf S_{3}$ remains valid. Note that the implied coupling of the plastic mechanisms would have to be postulated in the phenomenological plasticity theory.

The picture is simpler in the case of a triangular lattice where the loading (2) from $\bf T_{1}$ to $\bf T_{2}$ produces a mixture of only two ’pure’ states $\bf T_{1}$, $\bf T_{2}$, see Fig.5(d). The latter can be interpreted as the activation of a single plastic mechanism, the one favored by the loading, see Movie S2 in com .

To conclude, our model shows that crystal plasticity naturally arises from nonlinear elasticity, if the tensorial symmetry of the crystal lattice is properly accounted for. The memory of the atomic lattice in such infinitely periodic Landau theory is present in the form of the information about the affine mappings that leave the energy density invariant. Athermal evolution in the regularized model of this type can lead to temporal and spatial complexity which our analysis predicts to be highly sensitive to both, the crystallographic symmetry and the orientation of the crystal Weiss et al. (2015); Sparks and Maaß (2018). In particular, our study highlights the crucial role played in plastic deformation by the degenerate saddle points of the energy representing seemingly irrelevant, unstable crystallographic phases; for similar effects in other fields see Rehbein and Carpenter (2011); Wales (2003); Shtyk et al. (2017); Bociort and van Turnhout (2005). More generally, the proposed Landau theory perspective on crystal plasticity promises to become an important new tool in the study of inelasticity at the micro/nano scales where the conventional engineering theories fail to access strength, account for fluctuations and adequately describe size effects Greer and De Hosson (2011); Zhang et al. (2016, 2017).