Parallel three-dimensional simulations of quasi-static elastoplastic solids
PParallel three-dimensional simulations of quasi-static elastoplasticsolids
Nicholas M. Boffi a , Chris H. Rycroft a,b a John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02139 b Computational Research Division, Lawrence Berkeley Laboratory, Berkeley, CA 94720
Abstract
We extend to three dimensions a new projection method for simulating hypo-elastoplastic solidsin the quasi-static limit. The method is based on a surprising mathematical correspondenceto the incompressible Navier–Stokes equations, where the projection method of Chorin (1968)is an established numerical technique. We first review the method and present its extensionto three dimensions. We discuss the development of a three-dimensional parallel geometricmultigrid solver employed to solve a linear system for the quasi-static projection. We test themethod by simulating three-dimensional shear band nucleation and growth, a precursor tofailure in many materials. Although the method can be applied to any elastoplasticity model,we employ a physical model of a bulk metallic glass based on the shear transformation zonetheory. We consider several examples of three-dimensional shear banding, and examine shearband formation in physically realistic materials with heterogeneous initial conditions.
Keywords: elastoplasticity, Chorin-type projection method, multigrid methods, parallelcomputing, strain localization
1. Introduction
Modeling failure of materials is a fundamental problem in modern engineering and science.At small loads, the strain in a material is typically smooth, but as the loading is increased,the strain may become highly localized, which ultimately leads to failure [1]. This localizationof strain may be driven by geometrical effects such as necking [2, 3, 4, 5], or by instabilitiesin the material response, such as due to plastic yielding. In the latter case, positive feedbackcauses regions that have already plastically deformed to yield further, often resulting in alocalized plane of strain called a shear band [6, 7, 8, 9, 10]. Shear bands have been studiedanalytically [11, 12], but it remains challenging to model their formation in inhomogeneous,three-dimensional material samples, and numerical tools must be employed.Numerically modeling materials can be performed using the framework of elastoplasticity.For stress levels below the material yield stress, an elastoplastic material deforms purely
Email addresses: [email protected] (Nicholas M. Boffi), [email protected] (Chris H.Rycroft) a r X i v : . [ phy s i c s . c o m p - ph ] A ug lastically and returns to its undeformed state upon removal of the load; when the yieldstress is reached, the material begins to deform plastically, leading to permanent, irreversibledeformation that persists beyond load removal [13]. Elastoplasticity is complex and it admitsa number of mathematical descriptions [14, 15], each of which amounts to a specification of theinteraction of the elastic and plastic components of deformation at a microscopic level [16, 17].In this paper, we explore the hypo-elastoplasticity model [18], in which the Eulerian rate ofdeformation tensor is decomposed additively into elastic and plastic parts, D = D el + D pl [19].Hypo-elastoplasticity has some drawbacks, but it is well-suited to problems with small elasticdeformation and large plastic deformation. This makes it particularly appropriate for studyingfailure and shear banding in hard materials.The hypo-elastoplastic formulation has several numerical advantages. Since it is based onthe Eulerian rate of deformation tensor, it is well-suited to a fixed-grid framework. Fixedgrids have simpler topologies than their Lagrangian counterparts, and are easier to programand parallelize. This is particularly important in three dimensions, where the computationalexpense mandates parallelization, and where the implementation difficulty increases relativeto lower dimensional simulations. Fixed-grid methods are also the methods of choice for fluidsimulation [20, 21, 22], and they allow a wider range of numerical linear algebra techniquesto be used, such as the geometric multigrid method [23, 24].The additive decomposition of D , coupled with the linear-elastic constitutive relationand a continuum formulation of Newton’s second law, leads to a closed hyperbolic system ofpartial differential equations for the material velocity, stress, and variables intrinsic to theplasticity model. To properly resolve elastic waves, the timestep ∆ t for an explicit numericalscheme is restricted by the well-known Courant–Friedrichs–Lewy (CFL) condition [25]. TheCFL condition states that ∆ t must obey a constraint for numerical stability: ∆ t ≤ hc e , where c e is a typical elastic wave speed in the medium and h is the grid spacing.In metals and other materials of interest, elastic waves can travel at kilometers per second,while in many practical scenarios the timescale of loading is much longer than the elasticwave travel time [26]. The CFL condition thus poses a prohibitive limit on the timestep forprobing realistic timescales and system sizes, and the development of alternative simulationapproaches which avoid resolving elastic waves is necessary. By looking at the limit of longtimes and small velocities, one can show that Newton’s second law can be replaced by aconstraint that the stresses must remain in quasi-static equilibrium. However, these equationsare no longer a hyperbolic system, since the ability to time-integrate the velocity field is lost.Recently, Rycroft et al. [15] demonstrated a surprising mathematical analogy betweenquasi-static elastoplasticity and the incompressible Navier–Stokes equations for fluid flow.The incompressible Navier–Stokes equations combine an explicit equation for the velocity –dependent on the fluid pressure – with a requirement that the velocity must be divergencefree. Much like in the case of quasi-static elastoplasticity, the divergence-free requirementon the velocity field is obtained as a limit of an explicit equation for the pressure. A well-known algorithm for this setting is the projection method of Chorin [27, 28]. In Chorin’smethod, the velocity field is first updated explicitly, but this intermediate velocity does notobey the incompressibility constraint. An elliptic problem is solved that simultaneously2nforces the incompressibility constraint and enables computation of the pressure. Rycroft et al. translated Chorin’s projection method to quasi-static elastoplasticity by exploitingthe mathematical analogy between the two sets of limiting equations. The new methodwas quantitatively tested on some simple two-dimensional systems using a multi-threadedsingle-processor implementation. It was shown to be practical and efficient, enabling the useof timesteps that are orders of magnitude larger than those required by the CFL condition.In this paper, we extend the quasi-static projection method to three dimensions. Due tothe greatly increased problem size, we develop a distributed parallel implementation, whichcreates algorithmic challenges. In particular, the projection method requires that we solvea coupled elliptic equation for the components of velocity. We develop a custom parallelgeometric multigrid code to solve the resulting linear system.As a physical testbed for our methodology, we employ an athermal formulation of theshear transformation zone (STZ) theory of Falk, Langer, Bouchbinder and coworkers as aplasticity model for a bulk metallic glass [29, 30, 31, 32]. Metallic glasses naturally lendthemselves to study through the hypo-elastoplasticity framework, as their elastic deformationis generally small and well-described by a linear theory, yet they can exhibit significant plasticdeformation [33]. Their elastic moduli are typically on the order of 10–100 GPa, and henceexperimental loading conditions often place samples in the quasi-static regime [34]. They alsopresent interesting and poorly understood fundamental physics [35, 36, 37, 38], which can bedifficult to probe in dimensions greater than one without the numerical methods presentedhere.Simulation studies of bulk metallic glasses have been essential to our understanding oftheir properties. The development of the original STZ theory was guided by observationsof molecular dynamics simulations [39]. Simulations of the necking instability in a barunder uniaxial tension using the STZ plasticity model highlight the interplay between elasticand plastic deformation [40, 41, 42]. The physical mechanisms of fracture in BMGs wereexplored using the two-dimensional projection method [43], subsequently allowing the fracturetoughness of BMGs to be predicted across a wide range of experimental conditions [44].Later experimental measurements due to Ketkaew et al. demonstrated that these simulation-based predictions were quantitatively correct [45]. Indeed, testable predictions for complexamorphous systems such as BMGs are rare, and the development of efficient numericalmethods such as the ones presented here provide a way to generate them and to guide futureexperimental inquiry.The STZ theory is a useful test case for our method both physically and numerically, butthe numerical methodology is general and can be used for many plasticity models within thehypo-elastoplasticity framework. These could include free-volume based models of BMGs [46],plasticity models based on the random first-order transition theory of the glass transition [47],hypo-elastic materials [48, 49, 50], geophysical models [51, 52], and rate-independent plasticitymodels [53, 54, 55, 56]. We also emphasize that Chorin’s projection method represents a firststep towards more complex projection-based algorithms such as gauge methods [57, 58, 59] andpressure-Poisson methods [60, 61], and that we have laid the groundwork here to generalizethese algorithms to the case of hypo-elastoplasticity.3nder loading, BMGs exhibit shear bands [62, 63], which rapidly lead to material failure [26,64, 65] and are one of the primary limitations in employing BMGs in applications [66].Analytical work probing shear bands in amorphous materials is difficult, particularly in two orthree dimensions, which highlights a need for computational investigations. The developmentof our method enables the study of shear-banding in three dimensions at large scale andhigh resolution without excessive computational expense. Our simulations demonstrate thatthis scale and resolution is indeed necessary, and expose interesting fine-scale and uniquelythree-dimensional features of shear banding. Our methodology opens the door to futurestudies probing the shape, structure, and topology of shear bands, as well as the mechanismand statistical properties of their formation.The structure of this paper is as follows. In Section 2, we describe the relation betweenChorin’s projection method and the projection algorithm for hypo-elastoplasticity employedhere. In Section 3, we describe a finite-difference implementation of our projection method,and describe a forward-Euler based explicit method for solving the hypo-elastoplastic equationsin the non-quasi-static limit. We also discuss our parallel multigrid implementation used forthe stress projection. In Section 4, we demonstrate convergence between the explicit andprojection methods in a regime in which the two are expected to produce similar results, andstudy several interesting examples of shear banding dynamics in a metallic glass.
2. Projection methods for fluid dynamics and hypo-elastoplasticity
We consider an elastoplastic material with Cauchy stress tensor σ ( x , t ) and velocity field v ( x , t ) at a position x and time t . The total rate of deformation tensor D is defined as thesymmetric part of the velocity gradient, D = ( L + L T ) with L = ∇ v . For any field f ( x , t ),we define the advective time derivative by dfdt = ∂f∂t + ( v · ∇ ) f . The fundamental assumptionof hypo-elastoplasticity is that the rate of deformation tensor can be additively decomposedinto a sum of elastic and plastic parts, D = D el + D pl .For stiff elastoplastic materials with small elastic deformation, the linear elastic constitutivelaw provides an accurate description, D σ D t = C : D el = C : (cid:0) D − D pl (cid:1) . (1) C is the fourth-rank stiffness tensor, taken to be homogeneous and isotropic. With Lam´e’sfirst parameter λ and shear modulus µ , the components of C are given by C ijkl = λδ ij δ kl + µ ( δ ik δ jl + δ il δ jk ) [67]. The time derivative D σ D t = d σ dt − L σ − σ L T + Tr( L ) σ denotes theTruesdell objective stress rate.From Newton’s second law, the material velocity obeys the equation ρ d v dt = ∇ · σ (2)where ρ denotes the material density. Equations 1 & 2 form a hyperbolic system of equationsfor the stress and velocity fields, which can be solved explicitly using standard finite-difference4imulation methods. This hyperbolic system will resolve elastic waves, and so the timestep∆ t and grid spacing ∆ x must satisfy the CFL condition ∆ t ≤ ∆ x/c e for numerical stability,where c e is an elastic wave speed. In materials such as metals and metallic glasses, elasticwaves travel on the order of kilometers per second. Spatial discretizations capable of resolvingfine-scale features of interesting physical phenomena in these materials can be as small asmicrometers. For ∆ x = 1 µ m and c e = 1 km/s, the CFL condition requires ∆ t ≤ We consider a scenario in which plastic deformation occurs on a timescale much greaterthan the time for waves to propagate through the material. In this setting, macroscopicplastic deformation takes place due to the accumulation of small velocity gradients over longtimes. The details of a limiting procedure describing this physical regime were performed inprevious two-dimensional work [15] and will not be reproduced here.In this quasi-static limit, the equation for the velocity in Eq. 2 can be approximatelyreplaced by a constraint on the stress ∇ · σ ≈ . (3)Equation 3 dictates that each infinitesimal material element remains in approximate quasi-static equilibrium, and is thus referred to as the quasi-static constraint. The evolutionequation for the stress in Eq. 1 is unaffected by the limiting procedure, and hence Eq. 1 mustbe solved subject to the global constraint Eq. 3 to obtain solutions valid in this limit.At this stage, it is unclear how to do so. The velocity v appears in Eq. 1 through D , butthere is no longer an equation that can be integrated explicitly to solve for it. It is also notguaranteed that solutions of Eq. 1 subject to the constraint in Eq. 3 will agree with solutionsof Eq. 1 and Eq. 2. We now demonstrate an analogy between the computational issues presented in theprevious section and those encountered in incompressible fluid dynamics. Consider a fluidwith velocity v , pressure p , and density ρ . The fluid velocity field obeys the Navier–Stokesequation, d v dt = −∇ p + ν ∇ v . (4)While the fluid density satisfies dρdt = − ρ ( ∇ · v ) , (5)along with an equation of state linking the fluid density to the fluid pressure. Using anexplicit scheme to solve the hyperbolic system in Eqs. 4 & 5 will resolve sound waves in thefluid, which leads to timestep restrictions from the CFL condition. In the long-time limit,Eq. 5 is traded for the incompressibility constraint on the velocity field, ∇ · v = 0 . (6)5 igure 1: The projection-based timestepping scheme for (a) the velocity field in incompressible fluid dynamicsand (b) the stress tensor in quasi-static hypo-elastoplasticity. In both cases, an intermediate field value(denoted with a superscript ∗ ) is first computed which does not obey the divergence-free constraint. Thisintermediate field value is then projected back onto the manifold of divergence-free solutions to compute thefield at the next timestep. This limit reduces the coupled partial differential equations for the pressure and velocityto a single constrained equation for the velocity. The pressure is present in the equationfor the velocity, though its evolution equation has been exchanged for the incompressibilityconstraint; this is much like the quasi-static limit of hypo-elastoplasticity described in thepreceding section.Chorin [27, 28] developed a numerical method for this system of equations that involves theuse of an orthogonal projection, spurring significant research into related algorithms [60, 61].Such projection methods proceed via a two-step procedure, where an intermediate velocity v ∗ is first computed which does not obey the incompressibility constraint. v ∗ is then orthogonallyprojected onto the manifold of divergence-free solutions through the solution of an ellipticproblem for an auxiliary field related to the pressure. The process of projection simultaneouslyenforces the constraint and enables computation of the pressure field.One typical approach is to employ a Hodge decomposition [60, 61], v ∗ = v + ∇ φ, (7)where v is the desired divergence-free velocity field and φ is an auxiliary field. One thenupdates v ∗ for a fixed interval of time via the equation v ∗ t + ( v · ∇ ) v + ∇ q = ν ∇ v ∗ , (8)where ∇ q is an approximation to the pressure gradient. Substituting Eq. 7 into Eq. 8 leadsto a formula for the pressure ∇ p = ∇ ( q + φ t ) − ν ∇ ∇ φ, (9)from which p can be computed. The divergence of Eq. 7 implies that ∇ · v = 0 if φ is suchthat ∇ φ = ∇ · v ∗ . (10)6he Poisson problem in Eq. 10 can be solved for φ using standard techniques of numericallinear algebra, and the projection can be completed by computing v = v ∗ − ∇ φ. (11)Boundary conditions on φ in Eq. 10 depend on the physical scenario of interest, and are criticalfor obtaining higher-order methods [60, 61]. The algorithm typically proceeds by advancingEq. 8 for a single timestep, computing v according to Eq. 11, and then setting v ∗ = v [68].This procedure is schematically represented in discrete-time in Fig. 1(a). Projection methodsavoid the CFL condition associated with compressive waves in the fluid, and hence can usesignificantly larger timesteps than explicit methods. An extension of projection methodsknown as gauge methods do not reset v ∗ = v at the end of timestep, and instead allow it tocontinue to evolve during the computation [57, 58, 69, 70, 71]It is possible to demonstrate that solving Eq. 10 represents an orthogonal projection. Wedefine the inner product between two vector-valued fields, (cid:104) v , u (cid:105) = (cid:90) Ω v · u d x , (12)where Ω is the simulation domain. Using this inner product, we can compute (cid:104) v n +1 − v n , v n +1 − v ∗ (cid:105) = − (cid:90) Ω (cid:0) v n +1 − v n (cid:1) · ∇ φ ( x ) d x = (cid:90) Ω (cid:0) ∇ · v n +1 − ∇ · v n (cid:1) φ ( x ) d x = 0 , (13)thereby establishing that the projection v n +1 − v ∗ is orthogonal to the difference betweenthe two velocity fields, v n +1 − v n . We now formulate a three-dimensional projection method for solving Eq. 1 subject to thequasi-static constraint Eq. 3. We define an intermediate stress σ ∗ = σ + C : ∇ Φ , (14)where Φ ( x , t ) is an auxiliary vector field. We can solve for σ ∗ by dropping the C : D term inEq. 1, σ ∗ t + ( v · ∇ ) σ = L σ + σ L T − Tr( L ) σ + C : (cid:0) ∇ q − D pl (cid:1) . (15)In Eq. 15, q represents an approximation to the velocity v . Substituting Eq. 14 into Eq. 15,we find C : D = C : ∇ ( q − Φ t ) , (16)from which D can be computed. Taking the divergence of Eq. 14 and requiring ∇ · σ = , Φ must satisfy the equation ∇ · ( C : ∇ Φ ) = ∇ · σ ∗ . (17)7quation 17 is a linear system with source term ∇ · σ ∗ that can be solved for Φ . Hence,we can devise a projection method for quasi-static hypo-elastoplasticity by evolving Eq. 15for a single timestep, finding Φ according to Eq. 17, and projecting σ ∗ onto the manifold ofdivergence-free solutions by computing σ = σ ∗ − C : ∇ Φ . The algorithm may then proceedby setting σ ∗ = σ . We represent this algorithm schematically in Fig. 1(b). A projectionmethod is defined by the choice of the approximate velocity field q , the auxiliary vector field Φ , and the integration method for Eq. 15. By instead allowing σ ∗ to evolve over the course ofthe computation rather than setting σ ∗ = σ after each timestep, we expect it to be possibleto develop gauge methods for quasi-static hypo-elastoplasticity, in an analogous manner tothe case of fluid dynamics.As in the case of fluid dynamics, we can show that the projection is orthogonal in asuitable inner product. To do so, we define an inner product between two stress tensors asin [15], (cid:104) σ , σ (cid:48) (cid:105) = (cid:90) Ω σ : S : σ (cid:48) d x , (18)where S = C − is the stiffness tensor. Equation 18 computes the elastic strain energy of amaterial with stress field σ and strain field S : σ (cid:48) , or vice-versa by symmetry. Because S is asymmetric positive definite tensor for physically realistic Lam´e parameters, this definition isan inner product. By explicit computation, (cid:104) σ n +1 − σ n , σ n +1 − σ ∗ (cid:105) = (cid:90) (cid:0) σ n +1 − σ n (cid:1) : S : C : ∇ Φ d x = (cid:90) (cid:0) σ n +1 − σ n (cid:1) : ∇ Φ d x = − (cid:90) (cid:0) ∇ · σ n +1 − ∇ · σ n (cid:1) · Φ d x = 0 . (19)
3. Numerical implementation
In this section, we describe an implementation of an explicit forward Euler method tosolve Eqs. 1 & 2, as well as a specific instance of the quasi-static projection method in Eqs. 15& 17. We model elastoplastic deformation in a BMG using an athermal variant of the STZtheory.
As a plasticity model for a metallic glass, we use an athermal form of the STZ theorysuitable for studying diverse materials including BMGs below the glass transition temperature,dense granular materials, and soft materials such as foams or colloidal glasses [30, 31]. Withinthe STZ theory, irreversible molecular rearrangements are assumed to occur sporadicallythroughout an otherwise elastic material, and each rearrangement induces a small incrementof strain. The accumulation of many such events leads to macroscopic plastic deformation.These rearrangements are assumed to occur at rare, localized, sites known as STZs when localstresses surpass the material yield stress s y . Thermal fluctuations of the atomic configuration8 able 1: Material parameters used in this study, for both linear elasticity and the STZ model of amorphousplasticity. The Boltzmann constant k B is used to convert energetic values to temperatures. Parameter ValueYoung’s modulus E
101 GPaPoisson ratio ν K
122 GPaShear modulus µ ρ − Shear wave speed c s .
47 km s − Yield stress s Y τ − sTypical local strain (cid:15) c /k B Thermodynamic bath temperature T
400 KSteady state effective temperature χ ∞
900 KSTZ formation energy e z /k B configurational subsystem governing the rearrangements that occur at the STZs, and a kinetic/vibrational subsystem governing the thermal vibrations of atoms in their cage ofnearest neighbors [72].STZs may be conceptualized as clusters of atoms predisposed to configurational rearrange-ments when subjected to external shear [30]. Each rearrangement corresponds to a transitionin the configurational energy landscape; these transitions are usually towards a lower-energyconfiguration, but there is a small probability for a reverse transition. Before the applicationof external shear, the material sample is at a local minimum. External shear alters the shapeof the energy landscape, and can make transitions to other states considerably more likely.The density of STZs in space follows a Boltzmann distribution in an effective disordertemperature denoted by χ [73, 74, 75, 76]. χ governs the out-of-equilibrium configurationaldegrees of freedom of the material and has many properties of the usual temperature: it ismeasured in Kelvin, and can be obtained as the derivative of a configurational energy withrespect to a configurational entropy [39]. χ is distinct from the thermodynamic temperature T , though it plays the same role for the configurational subsystem as T does for thekinetic/vibrational subsystem.The plastic rate of deformation tensor is proportional to the deviatoric part of the stresstensor σ = σ − I tr( σ ), so that D pl = D pl σ ¯ s . ¯ s is a local stress measure given by theFrobenius norm of the deviatoric stress tensor, ¯ s = (cid:80) ij σ ,ij . The magnitude of the plastic9ate of deformation is given by τ D pl = e − e z /k B χ C (¯ s, T ) (cid:16) − s Y ¯ s (cid:17) , (20)where τ is a molecular vibration timescale, e z is a typical STZ formation energy, and k B isthe Boltzmann constant. C (¯ s, T ) represents the total STZ transition rate. With R ( ± ¯ s, T )denoting the forward and reverse rates between two configurational states, the total transitionrate is C (¯ s, T ) = ( R (¯ s, T ) + R ( − ¯ s, T )). The transitions follow a linearly stress-biasedthermal activation process, R ( ± ¯ s, T ) = exp (cid:18) − ∆ ∓ Ω (cid:15) ¯ sk B T (cid:19) . (21)∆ is a typical energetic barrier for a transition, Ω is a typical STZ volume, and (cid:15) is atypical local strain due to an STZ transition. While thermal fluctuations are neglected in theathermal model, the thermodynamic temperature still sets the magnitude of transition ratesin the system. Using the form Eq. 21 yields the overall transition rate C (¯ s, T ) = e − ∆ /k B T cosh (cid:18) Ω (cid:15) ¯ sk B T (cid:19) . (22)The effective temperature satisfies [30, 77, 35, 36] c dχdt = (cid:0) D pl : σ (cid:1) s Y ( χ ∞ − χ ) + l ∇ · (cid:0) D pl ∇ χ (cid:1) . (23)Equation 23 consists of a term causing growth to an asymptotic value χ ∞ and a diffusiveterm with diffusion length scale l . Both saturation to χ ∞ and diffusion occur in response toplastic deformation. The term D pl : σ is the rate of energy dissipated by externally appliedmechanical work, so that STZs are created and annihilated proportional to this rate, and c is an effective heat capacity. Eq. 23 is thus essentially a heat equation, representing thefirst law of thermodynamics for the configurational subsystem [30]. The interdependence ofEqs. 20 & 23 enables the development of shear bands via positive feedback, as increasing χ also increases D pl [35, 36]. For the explicit method, we discretize Eqs. 1, 2, and 23 with a forward Euler step. Thisleads to the coupled set of discrete-time equations ρ v n +1 − v n ∆ t = − ( v n · ∇ ) v n + ∇ · σ n + κ ∇ v n , (24) σ n +1 − σ n ∆ t = − ( v n · ∇ ) σ n + L n σ n + σ n (cid:0) L T (cid:1) n + Tr( L n ) σ n + C : (cid:18) D n − D pl ¯ s n σ n (cid:19) , (25) c χ n +1 − χ n ∆ t = (cid:0) ( D pl ) n : σ n (cid:1) s Y ( χ ∞ − χ n ) + l ∇ · (cid:0) D pl ∇ χ n (cid:1) . (26)10he small viscous stress term κ ∇ v in Eq. 24 is artificially imposed for numerical stability [78],but is not needed in the quasi-static method. In three dimensions, this term induces arestriction on the timestep ∆ t ≤ h κ . Hence, if κ is viewed as a physical constant, thiscondition is more restrictive than the CFL condition. However, for stability, it is sufficientto choose κ as scaling linearly with the grid spacing, in which case the timestep restrictionscales in the same way as the CFL condition. We now formulate a specific three-dimensional projection method for solving Eq. 1 subjectto the quasi-static constraint Eq. 3. We first neglect the C : D term in Eq. 1 and computean intermediate stress σ ∗ , σ ∗ = σ n + ∆ t (cid:18) L n σ n + σ n (cid:0) L T (cid:1) n − Tr( L n ) σ n − C : (cid:18) D pl ¯ s n σ n (cid:19)(cid:19) . (27)If the velocity at the next timestep v n +1 were known, we could compute D n +1 = 12 (cid:16)(cid:0) ∇ v n +1 (cid:1) + (cid:0) ∇ v n +1 (cid:1) T (cid:17) (28)and complete the forward Euler step in Eq. 27 as σ n +1 = σ ∗ + ∆ t (cid:0) C : D n +1 (cid:1) . (29)Taking the divergence of Eq. 29 and rearranging terms leads to the equation∆ t ∇ · (cid:0) C : D n +1 (cid:1) = −∇ · σ ∗ . (30)Equation 30 is a linear system for the velocity v n +1 involving mixed spatial derivatives. Thesource term is given by −∇ · σ ∗ . After solving for v n +1 , it can be used to compute σ n +1 viaEq. 29. Through this process, σ ∗ is projected back to be divergence-free, arriving at σ n +1 .The mixed derivatives in Eq. 30 increase the complexity of the projection for hypo-elastoplasticity when compared to the Poisson problem in fluid dynamics, but Eq. 30 cannevertheless be solved rapidly via standard techniques of numerical linear algebra such as themultigrid method. The multigrid method relies on the Gauss–Seidel method for iterativesmoothing of the solution, and Gauss–Seidel smoothing is guaranteed to converge if eitherthe linear system is (A) weakly diagonally dominant, or (B) symmetric positive definite. Ingeneral the linear system in Eq. 30 will not satisfy condition A, but will satisfy condition B.Hence Gauss–Seidel smoothing is guaranteed to converge, which we use as a component in amultigrid method—details of this multigrid solver are presented later. A connection to thegeneral continuous-time framework presented in Sec. 2.4 is provided in Appendix A.11 .4. Discretization and finite difference stencils The evolution equation for the stress, Eq. 1, depends on spatial derivatives of the velocity,while the equation satisfied by the velocity, Eq. 2, depends on spatial derivatives of the stress.We exploit this structure through a staggered grid with uniform spacing ∆ x = ∆ y = ∆ z = h .The stress tensor σ and effective temperature χ are stored at cell centers and indexed byhalf-integers, while the velocity v is stored at cell corners and indexed by integers, as shownin Fig. 2, left.Let ( ∂f /∂x ) i,j,k denote the partial derivative of a field f with respect to x evaluated atgrid point ( i, j, k ). The staggered centered difference is (cid:18) ∂f∂x (cid:19) i + ,j + ,k + = 14 h (cid:16) f i +1 ,j,k − f i,j,k + f i +1 ,j +1 ,k − f i,j +1 ,k + f i +1 ,j,k +1 − f i,j,k +1 + f i +1 ,j +1 ,k +1 − f i,j +1 ,k +1 (cid:17) . (31)Equation 31 averages four edge-centered centered differences surrounding the cell centerand has a discretization error of size O ( h ). The derivative at a cell corner is obtained bythe replacement ( i, j, k ) → ( i − , j − , k − ). The diffusive term appearing in the velocityupdate in Eq. 24 is computed via the standard centered difference formula, (cid:18) ∂ f∂x (cid:19) i,j,k = f i +1 ,j,k − f i,j,k + f i − ,j,k h . (32)The advective derivatives in Eqs. 24 & 25 must be upwinded for stability; we use the second-order essentially non-oscillatory (ENO) scheme [79]. With [ f xx ] i,j,k denoting the secondderivative with respect to x of the field f at grid point ( i, j, k ) computed using Eq. 32, theENO derivative is defined in the x direction as (cid:18) ∂f∂x (cid:19) i,j,k = 12 h − f i +2 ,j,k + 4 f i +1 ,j,k − f i,j,k if u i,j,k < (cid:12)(cid:12)(cid:12) [ f xx ] i,j,k (cid:12)(cid:12)(cid:12) > (cid:12)(cid:12)(cid:12) [ f xx ] i +1 ,j,k (cid:12)(cid:12)(cid:12) ,3 f i,j,k − f i − ,j,k + f i − ,j,k if u i,j,k > (cid:12)(cid:12)(cid:12) [ f xx ] i,j,k (cid:12)(cid:12)(cid:12) > (cid:12)(cid:12)(cid:12) [ f xx ] i − ,j,k (cid:12)(cid:12)(cid:12) , f i +1 ,j,k − f i − ,j,k otherwise . (33)Equation 33 uses the curvature of f to switch between an upwinded three-point derivativeand a centered difference. Versions of Eqs. 31, 32, & 33 in the y and z coordinates areobtained analogously.To solve Eq. 30 via numerical linear algebra, the spatial derivatives must first be discretizedusing finite differences. In addition to the finite differences discussed above, Eq. 30 alsocontains mixed partial derivatives. The xy -derivative is computed numerically as (cid:18) ∂ f∂x∂y (cid:19) i,j,k = f i +1 ,j +1 ,k − f i +1 ,j − ,k − f i − ,j +1 ,k + f i − ,j − ,k h , (34)with analogous expressions for other mixed partials.12 .5. Parallelization via MPI and domain decomposition We solve Eqs. 24 & 25 in parallel with a custom C++ code that uses the MPI libraryfor parallelization [80]. We use a global grid comprised of Q × M × N grid cells. Eachgrid cell is labeled by its lower corner index ( i, j, k ) as shown in Fig. 2, left, with indices inthe ranges i ∈ { , , . . . , Q − } , j ∈ { , , . . . , M − } , k ∈ { , , . . . , N − } . The globalgrid is split across processors using a standard Cartesian decomposition into P x × P y × P z subdomains. Specifically, each processor is indexed as ( I, J, K ) with I ∈ { , . . . , P x − } , J ∈ { , . . . , P y − } , and K ∈ { , . . . , P z − } . Processor ( I, J, K ) is responsible for the gridcells i ∈ {(cid:98) IQ/P x (cid:99) , . . . , (cid:98) ( I + 1) Q/P x (cid:99) − } ,j ∈ {(cid:98) J M/P y (cid:99) , . . . , (cid:98) ( J + 1) M/P y (cid:99) − } ,k ∈ {(cid:98) KN/P z (cid:99) , . . . , (cid:98) ( K + 1) N/P z (cid:99) − } . (35)Hence the subdomain assigned to each processor is identical up to one grid point in eachdirection.The finite difference stencils in Eqs. 31 & 32 require data from adjacent gridpoints, andthe ENO derivative in Eq. 33 requires data from at most two grid points away. On gridpoints within two points of a subdomain boundary, the derivative calculation can thereforerequire inaccessible data in a distributed memory setting. To handle this, we pad eachprocessor subdomain with ghost regions of width two. A ghost region is a cubical shell ofnon-physical grid points whose field values are filled with data from adjacent processors,so that each subdomain can freely and locally access information computed and storedin adjacent subdomains (Fig. 2, right). At the simulation boundaries, ghost regions areused to enforce boundary conditions. For periodic boundary conditions, the ghost regionswrap around to processor subdomains on the other side of the simulation. For non-periodicconditions, the ghost grid points store values which enforce the desired boundary conditions.At the start of each timestep, each processor communicates with 26 nearby processors vianon-blocking communication, sharing six faces, twelve edges, and eight corner regions. Eachprocessor sends data to nearby processors, receives the data it requires from the same nearbyprocessors, and loads that data into its ghost regions. The total cost of parallel communicationscales with the surface area of a processor subdomain. Usually, the total number of processors P tot is chosen to be power of two, and the precise subdomain decomposition is determinedby considering all triplets ( P x , P y , P z ) that satisfy P x P y P z = P tot , and selecting the one thatminimizes the surface area between the subdomains. We solve Eq. 30 for the velocity using a custom parallel implementation of the geometricmultigrid method, a multi-resolution linear system solver that is particularly suited to ellipticproblems that take place on a physical grid [23]. Let G be the original grid, and let A x = b be the linear system to solve on this grid. In the multigrid method, a hierarchy ofprogressively coarser grids G , G , . . . , G g is introduced. In our implementation, if G k hasresolution Q k × M k × N k , then G k +1 has resolution (cid:100) Q k / (cid:101) , (cid:100) M k / (cid:101) , (cid:100) N k / (cid:101) . Interpolation13 igure 2: (Left) The arrangement of discretized simulation fields. The cube here corresponds to a single gridcell of side length h . Velocities v are stored at cell corners denoted by black spheres and indexed by integers.Stresses σ and effective temperatures χ are stored at cell centers indicated by the purple sphere and indexedby half integers. (Right) Processor ghost regions. The two solid cubes represent two adjacent processorsubdomains; the boundary between them is indicated by a gray plane in the center of the figure. Surroundingeach processor subdomain is a transparent two-grid-point ghost region bounded by black dashed lines. Forclarity, the ghost region grid in the left processor subdomain has been drawn in a thin gray, while the ghostregion grid in the right processor subdomain has been drawn in black. These two processors communicatethe overlapping rectangular strip surrounding the separating plane in the center of the figure. operators T k : G k → G k − are introduced based on linear interpolation, and restrictionoperators R k : G k → G k +1 are introduced based on local averaging. Both T k and R k can berepresented as rectangular matrices, and in our implementation R k − = T T k —this conditionis not necessary for a practical implementation, but is useful in some convergence proofs [24].Our multigrid implementation uses the standard V-cycle [23, 24] with two pre-smoothingsteps and two post-smoothing steps. On G the grid is decomposed among the processors inthe same way as the simulation fields (Fig. 2(b)). The smoothing steps are performed usingthe Gauss–Seidel method on each processor, with the ghost regions being synchronized aftereach step. This requires building a representation of the linear system on each grid, whichwe do via recursive matrix multiplication [81, 82], A k = R k − A k − T k . (36)The implementation works with periodic and non-periodic boundary conditions, and arbitrarygrid dimensions. As the grids are coarsened, the amount of work on each grid is rapidlyreduced, to the point where it is no longer effective for all processors to share the work. Theimplementation therefore has the ability to amalgamate the coarser problem onto a smaller14et of processors, with the rest remaining idle.The multigrid implementation uses C++ templates, so that the linear system can becompiled to work with an arbitrary data type. For the current problem, b is given by thesource term − ∆ t ∇ · σ and x contains values of v n +1 across the entire grid. Hence, wecompile the multigrid library where the elements of b and x are 3-vectors, and the elementsof A are 3 × A is sparse, and a grid point ( i, j, k ) isonly coupled to the 27 grid points in the 3 × × i + {− , , } , j + {− , , } , k + {− , , } ) in our discretization scheme.
4. Shearing between two parallel plates
In the following sections, we consider several material samples being sheared between twoparallel plates. This example is experimentally relevant, has simple boundary conditions,demonstrates complex shear banding dynamics [65, 63, 62, 26, 66, 38, 37, 64], and has beenstudied previously in two dimensions [15]. It represents a natural physical scenario to comparethree-dimensional results to two-dimensional results, compare simulation data to experiments,and to quantitatively compare the explicit and quasi-static methods.The domain occupies − L ≤ x < L, − L ≤ y < L , and − γL ≤ z ≤ γL with γ = and L = 1 cm. A natural unit of time is given as t s = L/c s where c s = (cid:112) µ/ρ is the materialshear wave speed, and we measure time in this scale. In all simulations, we consider adomain periodic in the x and y directions with shear velocity applied on the top and bottomboundaries in z . The boundary conditions are given by v ( x, y, ± γL, t ) = ( ± U ( t ) , , , (37)where the function U ( t ) is given by U ( t ) = (cid:40) U B tt s if t < t s , U B otherwise. (38)The ramp-up in the function U ( t ) prevents a large deformation rate near the boundarythat would be present with U ( t ) = U B immediately at t = 0. The elasticity and plasticityparameters are defined in Table 1. From these values, the natural timescale is t s = 4 . µ s.A diagram of the global three-dimensional grid and the ghost regions at simulation bound-aries used for implementing the boundary conditions is shown in Fig. 3. The cell-cornered gridpoints run according to i ∈ { , . . . , Q − } , j ∈ { , . . . , M − } , and k ∈ { , . . . , N } ; becausethe grid is non-periodic in z there is an extra grid point in this direction. The velocities atgrid points with k = 0 and k = N are fixed according to the boundary velocity in Eq. 37.The cell-centered grid points run according to i ∈ { , , . . . Q − } , j ∈ { , , . . . M − } , and k ∈ { , , . . . N − } . Ghost layers of cell-centered grid points are at ( i, j, − ), ( i, j, − ),( i, j, N + ), and ( i, j, N + ). The values of σ and χ in the ghost layers are linearly inter-polated from the two most adjacent layers, to ensure that these fields remain free on theboundary. At the simulation boundaries in the x and y directions, ghost points outside thesimulation domain are filled with values that wrap around.15 igure 3: A diagram of the simulation grid layout for the simplified case of ( Q, M, N ) = (2 , , ± x and ± y faces and are surrounded bysee-through rectangular prisms. For clarity, these are omitted from the ± z faces. Corner-centered ghostpoints and cell-centered ghost points are smaller and are shown in lighter pink and green than their physicalcounterparts. The ghost points adjacent to the ± x and ± y faces wrap around, and are used to enforceperiodic boundary conditions. In the z direction, the ghost grid points are used to linearly interpolate the σ and χ values, leaving both fields free on the boundary. In the z direction, there is one extra corner-centeredgrid point, giving the appearance of a grid of size 2 × ×
3. This grid point is used to enforce shear boundaryconditions on the velocity field, but the equivalent cell-centered grid points are used to store ghost σ and χ values. In each case considered in the following sections, the physics of the material sample isencoded in the initial effective temperature distribution χ ( x , t = 0). Following the introductionin Sec. 3.1, χ is a continuum-scale variable that encodes the density of STZs, and hence itsinitial condition affects the future evolution of plastic deformation within the material. We now demonstrate the qualitative equivalence between results computed with theexplicit and quasi-static methods. We consider an initial condition corresponding to a finitecylindrical inclusion χ ( x , t = 0) = 600 K + (200 K) e − z + y ) /L , (39)for x > − L/ x < L/
2, and 600 K otherwise. Initially the cylindrical inclusion is slightlymore amenable to plastic deformation, and hence we expect to see a shear band nucleatefrom it. To visualize the effective temperature field in three dimensions, we use a custom16 igure 4: The effective temperature field χ in the initial configuration for the quasi-static to explicit methodquantitative comparison. Here, a = 0 .
25 and η = 0 . opacity function defined as O ( x ) = (cid:16) χ ( x ) − χ bg χ ∞ − χ bg (cid:17) if χ ( x ) − χ bg χ ∞ − χ bg > , exp (cid:16) − a (cid:16) χ ∞ − χ bg χ ( x ) − χ bg (cid:17) η (cid:17) otherwise. (40)Equation 40 sets the opacity of a grid point based on the value of χ ( x ). The parameters a and η are chosen on a case-by-case basis to reveal the most interesting features . The initialcondition is depicted in Fig. 4. The grid is of size 64 × ×
32 to accommodate the limitationsof the explicit simulation method, corresponding to a grid spacing h = L . A viscous stressconstant of κ = 0 . L t s is used in Eq. 24. The timestep for the explicit method is ∆ t e = h t s L .A typical applied shear velocity that is comparable to a realistic loading rate in alaboratory experiment is U b = 10 − L/t s [65]. With this velocity, running an explicit simulationis prohibitively expensive due to the CFL condition. To ensure that significant plasticdeformation occurs on timescales reachable by the explicit method, a scaling parameter ζ is introduced. The molecular vibration timescale τ is rescaled to τ ζ − and the appliedshear velocity is inversely rescaled to U B = 10 − L/t s ζ . The simulation is conducted until afinal time of t f = 2 × t s /ζ . As ζ approaches zero, the quasi-static limit of Eqs. 1 & 2 isformally approached. We therefore expect greater agreement for smaller values of ζ . Due tothe appearance of τ in Eq. 20, the introduction of ζ has the effect of linearly scaling themagnitude of plastic deformation by a factor of ζ . A quasi-static timestep of ∆ t qs = 200 t s /ζ is used.In Fig. 5, we show three snapshots of the effective temperature field with ζ = 10 from Ideally, we would like to use the same opacity parameters for all plots. However, due to significantvariations in the ranges of the χ fields and in their spatial structures, we found it was necessary to set theparameters on an individual basis. We note, however, that the color scale is absolute across all simulations. t = 50 t s , t = 75 t s , and t = 100 t s respectively. Theexplicit simulation is shown on the left and the quasi-static simulation is shown on the right.The results are qualitatively similar in all three snapshots. At t = 50 t s , a shear band beginsto emerge, nucleating outwards from the center of the simulation. A thin region of higher χ is visible in the center of the band. By t = 75 t s , the shear band has fully formed and spansthe system. At t = 100 t s , the band grows stronger and χ continues to increase.Figure 6 shows cross-sections in z for fixed x = 0 and y = 0 of (cid:107) σ (cid:107) qs − (cid:107) σ (cid:107) e for severaltime points before the onset of plastic deformation, highlighting some differences between thetwo methods. The explicit simulation exhibits oscillations due to elastic waves propagatingthrough the medium. Because the quasi-static method does not resolve these elastic waves,the oscillations are apparent in the deviatoric stress differences. When plastic deformationsets in, plasticity-induced damping removes the elastic waves and the agreement improves. Having demonstrated the qualitative agreement between the two simulation methodologiesfor ζ = 10 in the previous section, we now examine convergence as ζ is decreased. The samesimulation geometry, boundary conditions, and initial conditions in the effective temperaturefield are used here as in the previous section. To quantitatively compute the agreementbetween the explicit and quasi-static methods, we define a norm on simulation fields f , (cid:107) f (cid:107) ( t ) = (cid:115) γL (cid:90) γL − γL (cid:90) L − L (cid:90) L − L | f ( x , t ) | dx dy dz. (41)The integral in Eq. 41 runs over the entire simulation domain and is computed numericallyvia the trapezoid rule. The appearance of | · | in Eq. 41 is taken to be the Euclidean norm forvectors, absolute value for scalars, and the Frobenius norm for matrices.Equation 41 is evaluated for χ qs − χ e , σ qs − σ e , and v qs − v e at intervals of 0 . t s in pairsof explicit and quasi-static simulations with ζ = 10 , × , . × and 1 . × . Ineach case, the norm is non-dimensionalized using the quantities χ ∞ , s Y , and U B to ensure allvalues are of order unity. For each simulation, the explicit timestep is ∆ t e = t s h L and thequasi-static timestep is ∆ t qs = t s ζ .Plots of all three norm values are shown as a function of time in Fig. 7(a) for the valueof ζ = 1 . × . Oscillations due to elastic waves are visible in all simulation fields untilaround t ≈ t s when the yield stress is reached. After the onset of plastic deformation, thenorm in effective temperature increases steadily, most rapidly during the period of shearband nucleation from t ≈ t s to t ≈ t s . The disagreements in σ and v decrease duringthe elastic region, and steadily increase after plastic deformation begins.Figures 7(b), 7(c), and 7(d) show the quantitative comparisons as a function of time forvalues of ζ = 10 , × , . × , and 1 . × for v , χ , and σ respectively. In all plots,better agreement with smaller ζ is observed during the elastic regime and during the onset ofplasticity while t ≤ t s . After shear band nucleation from 12 t s ≤ t ≤ t s , all values of ζ have roughly equal error magnitudes, with slightly greater agreement for higher values of ζ .This is consistent with previous comparisons in two dimensions, where the dominant factor18 igure 5: Snapshots of the effective temperature distribution χ ( x , t ) for the explicit simulation (left) andquasi-static simulation (right) for ζ = 10 . The simulation fields are qualitatively similar. In all plots, a = 0 . η = 1 . t = 50 t s . (c,d) t = 75 t s . (e,f) t = 100 t s . igure 6: (a) The magnitude of the deviatoric stress tensor (cid:107) σ (cid:107) for the explicit and quasi-static simulationmethods along a cross section in z for x = 0 and y = 0 fixed. Results for the explicit and quasi-staticsimulation methods are shown in dashed and solid lines respectively. Oscillations at t = 5 t s and t = 10 t s aredue to elastic waves propagating through the medium in the explicit simulation, but are difficult to see byeye at this scale, see (b). As plasticity kicks in past t = 15 t s , these waves damp out. (b) The difference in themagnitude of the deviatoric stress tensor (cid:107) σ (cid:107) for the explicit and quasi-static simulation methods, along across section in z for x = 0 and y = 0 fixed. The oscillations are due to elastic waves propagating throughthe medium in the explicit simulation.Figure 7: L norm of the χ , v , and σ simulation field differences between the explicit and quasi-static methodcomputed using Eq. 41 and normalized by the respective characteristic variables. (a) A comparison of thefour different field norms, for the value of ζ = 10 . The remaining three panels show the differences in (b)velocity, (c) effective temperature, and (d) stress, respectively, for a range of values of ζ . ζ itself [15].A method to reduce the differences in discretization is to increase the background χ field.With higher values of background χ , finer-scale features in the shear banding dynamics areless prominent. This ensures that differences in the spatial discretization will be minimized.There is also less rapid development of the shear band, and thus the difference in timestepbetween the two methods will be less significant. Snapshots of the effective temperature fieldare shown in Fig. 8 at t = 10 t s for background χ field of χ bg = 600 K ,
650 K ,
700 K. Figure9 confirms that the differences between the two types of simulation decreases as χ bg increases. We now turn to simulating realistic physical timescales with the quasi-static method,where the scaling parameter is ζ = 1. We first consider the nucleation of shear bands fromlocalized imperfections of higher χ . Physically, this describes defects within the materialstructure which may be particularly susceptible to plastic deformation [83]. To begin, weconsider a single defect, corresponding to an initial χ field of the form χ ( x , t = 0) = 550 K + (170 K) exp (cid:18) − (cid:107) x (cid:107) L (cid:19) . (42)The simulation is performed on a grid of size 256 × × h = L/ l appearing in Eq. 23 is fixed at 3 h and sets the width ofthe shear bands. The boundary velocity is set to a value of U B = 5 × − L/t s , and thesimulation is conducted to a final time of t f = 4 × t s , using a quasi-static timestep of∆ t = 200 t s . For three-dimensional visualization, we use the opacity function from Eq. 40.Snapshots of the effective temperature field at various time points are shown in Fig. 10.The initial condition is shown in Fig. 10(a). At t = 10 t s in Fig. 10(b), the defect has startedto expand. By t = 2 × t s in Fig. 10(c), a shear band begins to nucleate, indicated by aquadrupolar structure emanating from the defect. The background χ field also begins toincrease, as demonstrated by the presence of the transparent light blue background. By t = 2 . × t s in Fig. 10(d), a distinct system-spanning band has become clear with apropagating front visible near its center. The band displays no curvature in either of the x or y directions. By t = 150 t s in Fig. 10(e), a prominent band has formed, and there is nolonger a visible propagating front. The band continues to grow stronger and thicker through t = 175 t s in Fig. 10(f) and t = 200 t s (not shown).We now introduce a second defect to highlight some three-dimensional characteristics ofshear banding. We expect that the relative size and displacement between the defects willdetermine the dynamics, with the possibility of forming a single shear band that connectsthe two. The initial effective temperature field is χ ( x , t ) = 550 K + (200 K) (cid:32) exp (cid:32) − (cid:107) ( x − X ) (cid:107) L (cid:33) + exp (cid:32) − (cid:107) ( x − X ) (cid:107) L (cid:33)(cid:33) . (43)Two cases of Eq. 43 are considered. First, we take X = ( − . , − . , .
35) and X =( − . , . , . y = 0 plane with the same21 igure 8: Qualitative demonstration of the effect of increasing the background χ field with χ bg set to (a/b)600 K (c/d) 650 K and (e/f) 700 K. All snapshots are displayed at t = 10 t s for a value of ζ = 10 , withfixed values of a = 0 .
75 and η = 3 in the opacity function. Simulation results are shown for the explicitmethod on the left and the quasi-static method on the right. For lower χ bg , the shear band is more prominent,develops more rapidly, and has more fine-scale features, ensuring that the differences in spatial and temporaldiscretizations become more pronounced. igure 9: Normalized L difference in (a) v (b) χ and (c) σ between the explicit and quasi-static methods forvarious choices of background χ field. Agreement improves as χ bg increases due to a reduction in fine-scalefeatures that differ between the two methods due to differences in the discretization. x coordinate, a slight offset in z , and different sizes. The results for this case using the samesimulation parameters as Fig. 10 are shown in Fig. 11. Second, we take X = ( − . , − . , . X = (0 . , − . , . x and y interchanged. The results for this case again with the same grid size, quasi-static timestep,and boundary velocity as for the single defect are shown in Fig. 12. The initial configurationsare shown in Figs. 11(a) and 12(a).There is an interesting contrast between the time sequences displayed in Figs. 11 and12. Much like in the single defect simulations, t = 50 t s in Figs. 11(b) and 12(b) displaysexpansion of the defects, and t = 100 t s in Figs. 11(c) and 12(c) shows the initiation of shearband nucleation. In Fig. 11(d), we see the formation of a single curved band connecting thetwo defects, while in Fig. 12(d), the band is flat. Figs. 11(e), 11(f), 12(e), and 12(f) makethis more clear as the band becomes more defined. The curvature seen in Fig. 11 is in thedirection orthogonal to shear.The dependence of band curvature on the relative orientation of the two defects can bebest understood in terms of the qualitative structure of Fig. 10(c). There is a substantialextension of elevated χ along the x direction (parallel to shear, in-plane), a small extensionalong the y direction (orthogonal to shear, in-plane), and a moderate extension along the z direction (orthogonal to shear, out-of-plane). In Fig. 11, the defects are offset in y and z . Because the χ field is stronger in z than in y , this relative placement of the defects canaccomodate curvature along the y direction. On the other hand, in Fig. 12, the defects areoffset in x and z . The strength of the χ field extension in the x direction is great enoughthat the flat, horizontal regions of the two forming bands reach each other. The two bandsjoin into one fatter flat band. 23 igure 10: Snapshots of the effective temperature distribution χ ( x , t ) for a quasi-static simulation with ζ = 1.The initial condition is given in Eq. 42, corresponding to a small Gaussian defect at the center of the material. a = 0 .
35 and η = 1 . a = 0 . η = 1 . t = 0 t s . (b) t = 10 t s . (c) t = 2 × t s . (d) t = 2 . × t s . (e) t = 3 × t s . (f) t = 4 × t s . igure 11: Snapshots of the effective temperature distribution χ ( x , t ) for a quasi-static simulation with ζ = 1.The initial condition is given in Eq. 43 with X = ( − . , − . , .
35) and X = ( − . , . , . a = 0 .
35 and η = 1 . a = 0 . η = 1 . t = 0. (b) t = 10 t s . (c) t = 2 × t s . (d) t = 2 . × t s . (e) t = 3 × t s . (f) t = 4 × t s . igure 12: Snapshots of the effective temperature distribution χ ( x , t ) for a quasi-static simulation with ζ = 1.The initial condition is given in Eq. 43 with X = ( − . , − . , .
35) and X = (0 . , − . , . a = 0 .
35 and η = 1 . a = 0 . η = 1 . t = 0. (b) t = 10 t s . (c) t = 2 × t s . (d) t = 2 . × t s . (e) t = 3 × t s . (f) t = 4 × t s . We now turn to a set of more complex initial conditions in the effective temperature field.Results for circular initial conditions parallel and perpendicular to the direction of shear areshown in Figs. 13 and 14 respectively. Mathematically, the initial conditions are d = (cid:114) y L + z L − ,χ ( x , t = 0) = 550 K + (200 K) exp (cid:0) − (cid:0) d + x (cid:1)(cid:1) , (44)and d = (cid:114) x L + z L − ,χ ( x , t = 0) = 550 K + (200 K) exp (cid:0) − (cid:0) d + y (cid:1)(cid:1) , (45)representing circles in the yz and xz planes respectively. Simulations were carried out usingthe same simulation geometry, discretization, quasi-static timestep, and boundary velocity asin the previous section. The initial conditions in Eqs. 44 and 45 are displayed in Figs. 13(a),14(a) respectively.By t = 8 × t s in Figs. 13(b) and 14(b), little has changed. At t = 1 . × t s inFigs. 13(c) and 14(c), differences due to the orientation of the circles become clear. The circleoriented perpendicular to shear closes vertically into a disk. The circle oriented along shearexhibits signatures of shear band nucleation at four equally space points. At t = 2 . × t s in Fig. 13(d), the disk has expanded and has developed a pointed front in the y direction.There are also two thick, curved bands forming close together near the center of the disk.In Fig. 14(d), there are two thinner, well-separated bands forming off the top and bottomof the circle with propagating fronts. By t = 3 . × t s in Figs. 13(e) and 14(e), thesedifferences have become even more prominent. The disk is still clear in Fig. 13(e) emergingfrom the two bands, and these bands are seen to have a curved structure in the y direction.They are also fatter and less separated than the bands seen in Fig. 14(e). These featurescontinue to develop into the final pane at t = 4 × t s . Taken together, Figs. 13 and 14demonstrate another example of the dependence of shear banding structure and dynamics onthe orientation of initial conditions in the χ field with respect to shear. In this section, we consider the case of a randomly distributed initial effective temperaturefield. The initial conditions presented in the previous sections provide insight into the27 igure 13: Snapshots of the effective temperature distribution χ ( x , t ) for a quasi-static simulation with ζ = 1.The initial condition is given by Eq. 44 with the circle oriented along the yz plane. a = 0 .
35 and η = 1 . a = 0 . η = 1 . t = 0 (b) t = 8 × t s (c) t = 1 . × t s (d) t = 2 . × t s (e) t = 3 . × t s (f) t = 4 × t s . igure 14: Snapshots of the effective temperature distribution χ ( x , t ) for a quasi-static simulation with ζ = 1.The initial condition is given by Eq. 45 with the circle oriented along the xz plane. a = 0 .
35 and η = 1 . a = 0 . η = 1 . t = 0 (b) t = 8 × t s (c) t = 1 . × t s (d) t = 2 . × t s (e) t = 3 . × t s (f) t = 4 × t s . χ is most faithful to this fundamental assumption [84]. Randominitial conditions are thus expected to shed the most light on the structure of shear bandsobserved in experiments. The randomly fluctuating χ field furthermore leads to the formationof multiple shear bands, and potentially enables the study of shear band interactions in theSTZ model [85].We first populate the grid and a shell of ghost points with random variables χ ζ ( x ) usingthe Box–Muller algorithm. With µ χ and σ χ respectively denoting the desired mean andstandard deviation, we perform the convolution χ ( x ) = σ χ N (cid:88) r ∈ V (cid:48) e − (cid:107) x − r (cid:107) l c χ ζ ( r ) + µ χ , N = (cid:115)(cid:88) r ∈ V e − (cid:107) r (cid:107) l c . (46)where V denotes the set of grid points and V (cid:48) denotes the set of grid points with the additionof the ghost points. Equation 46 ensures that the effective temperature value at each pointis normally distributed with mean µ χ and standard deviation σ χ . In practice, the sums inEq. 46 are performed with a cutoff length scale specified as a multiplicative factor of theconvolution length scale l c , and the number of ghost points in V (cid:48) is set by the choice of cutofflength scale. Results for a random initialization with µ χ = 550 K, σ χ = 15 K, l c = 10 h anda cutoff factor of 5 (leading to 50 ghost points in each direction for the convolution) areshown in Fig. 15. The grid is of size 768 × × t = 1 × t s in Fig. 15(b), the effective temperature has increased somewhat uniformlyacross the grid. At t = 4 × t s in Fig. 15(c), both horizontal and vertical shear bandsbegin to nucleate throughout the simulation. Slightly later at t = 6 × t s in Fig. 15(d), amultitude of thin, system-spanning horizontal bands connected by vertical bands have begunto emerge. Curvature is present in the horizontal bands both parallel and orthogonal to thedirection of shear. The shear bands become increasingly prominent and grow in number by t = 8 × t s and t = 10 t s in Figs. 15(e) and (f) respectively, where crossing and branchingpatterns in the many bands are observed.
5. Conclusion
Expanding on prior two-dimensional work [15], we have developed a three-dimensionalnumerical method for simulating quasi-static elastoplastic materials by analogy with theprojection method of Chorin for incompressible fluid dynamics [27, 28]. The method isparticularly suitable for stiff materials with small elastic deformation and high elastic wavespeeds, where plastic deformation is often quasi-static, and the hypo-elastoplastic assumptionis valid. In these materials, the timestep of an explicit method is limited by the CFL condition—making macroscopic timescales and realistic loading rates prohibitive to simulate—while thequasi-static method has no such restriction. 30 igure 15: Snapshots of the effective temperature distribution χ ( x , t ) for a quasi-static simulation with ζ = 1.(a) t = 0, a = 0 .
25, and η = 1 .
3. (b) t = 2 × t s , a = 0 .
25, and η = 1 .
3. (c) t = 4 × t s , a = 0 .
45, and η = 1 .
75. (d) t = 6 × t s , a = 0 .
5, and η = 1 .
8. (e) t = 8 × t s , a = 0 .
75, and η = 2 .
0. (f) t = 10 t s , a = 0 .
75, and η = 2 . χ , and high-resolution banding with a random initialization in χ .Essential to these studies was the development of a high performance, parallelized, three-dimensional geometric multigrid solver, along with an MPI and C++-based implementationof the projection algorithm. The method, though applied to the study of BMGs in this work,is independent of the plasticity model.The method is based on a surprising correspondence between the variables ( v , p ) inincompressible fluid dynamics and the variables ( σ , v ) in quasi-static hypo-elastoplasticity.As part of the development of the method, we introduced an auxiliary vector field Φ , whichplays a role analogous to the auxiliary scalar field φ used in projection algorithms for fluiddynamics. The choice of φ , along with careful consideration of its boundary conditions,leads to higher-order projection and gauge methods in fluid dynamics. The analogy betweenthe two auxiliary fields is reminiscent of the analogy between the physical fields that ledto the development of the quasi-static method, and it suggests that similar approachesmay generalize to quasi-static hypo-elastoplasticity [60, 57, 58, 59]. Gauge methods couplenaturally with discontinuous Galerkin discretization (dG) methods, and the theoreticaldevelopments presented here thus pave the way for the possibility of similar dG methods inquasi-static hypo-elastoplasticity, which may help resolve fine-scale features of instabilitiessuch as shear bands. The development of such higher-order methods is an interesting avenuefor future work. Appendix A. Connection to the continuous-time framework
We can make a connection to the general continuous-time framework presented in Sec. 2.4as follows. By comparison of Eqs. 17 and 30, we can identify Φ = ∆ t v . Equation 16 thensays that C : D n +1 = C : (cid:18) ∇ q + ∆ t ∂ D n +1 ∂t (cid:19) . (A.1)Recall that q is chosen to be the best available approximation to v n +1 , and note that bysymmetry of C , C : ∇ q = C : 12 (cid:16) ∇ q + ∇ ( q ) T (cid:17) . (A.2)Equation A.1 thus says the following: C : D n +1 is given by the best available guess beforethe solve for v n +1 - C : ∇ q - plus an O (∆ t ) correction constructed via a first-order Taylorexpansion in time. Acknowledgments
The authors thank Eran Bouchbinder and Avraham Moriel (Weizmann Institute of Science)for useful discussions about this work. This work was supported by the National ScienceFoundation under Grant Nos. DMR-1409560 and DMS-1753203. N. M. Boffi was supported32y a Department of Energy Computational Science Graduate Fellowship. C. H. Rycroftwas partially supported by the Applied Mathematics Program of the U.S. DOE Office ofAdvanced Scientific Computing Research under contract number DE-AC02-05CH11231.
References [1] M. Jir´asek, Objective modeling of strain localization, Revue Fran¸caise de G´enie Civil6 (6) (2002) 1119–1132. doi:10.1080/12795119.2002.9692735 .[2] J. Hutchinson, J. Miles, Bifurcation analysis of the onset of necking in an elastic/plasticcylinder under uniaxial tension, J. Mech. Phys. Solids 22 (1) (1974) 61–71. doi:https://doi.org/10.1016/0022-5096(74)90014-3 .[3] A. K. Ghosh, Tensile instability and necking in materials with strain hardening andstrain-rate hardening, Acta Metallurgica 25 (12) (1977) 1413–1424. doi:10.1016/0001-6160(77)90072-4 .[4] V. Tvergaard, Necking in tensile bars with rectangular cross-section, Computer Meth-ods in Applied Mechanics and Engineering 103 (1) (1993) 273–290. doi:10.1016/0045-7825(93)90049-4 .[5] H. Guo, P. F. Yan, Y. B. Wang, J. Tan, Z. F. Zhang, M. L. Sui, E. Ma, Tensileductility and necking of metallic glass, Nature Materials 6 (2007) 735–739. doi:10.1038/nmat1984 .[6] J. W. Rudnicki, J. R. Rice, Conditions for the localization of deformation in pressure-sensitive dilatant materials, Journal of the Mechanics and Physics of Solids 23 (6) (1975)371–394. doi:10.1016/0022-5096(75)90001-0 .[7] J. Hutchinson, V. Tvergaard, Shear band formation in plane strain, International Journalof Solids and Structures 17 (5) (1981) 451–470. doi:10.1016/0020-7683(81)90053-6 .[8] S. Harren, H. D`eve, R. Asaro, Shear band formation in plane strain compression, ActaMetallurgica 36 (9) (1988) 2435–2480. doi:10.1016/0001-6160(88)90193-9 .[9] R. Conner, Y. Li, W. Nix, W. Johnson, Shear band spacing under bending of zr-basedmetallic glass plates, Acta Materialia 52 (8) (2004) 2429–2434. doi:10.1016/j.actamat.2004.01.034 .[10] H. Bei, S. Xie, E. P. George, Softening caused by profuse shear banding in a bulk metallicglass, Phys. Rev. Lett. 96 (2006) 105503. doi:10.1103/PhysRevLett.96.105503 .[11] R. Asaro, J. Rice, Strain localization in ductile single crystals, J. Mech. Phys. Solids25 (5) (1977) 309–338. doi:10.1016/0022-5096(77)90001-1 .[12] P. Steif, F. Spaepen, J. Hutchinson, Strain localization in amorphous metals, ActaMetallurgica 30 (2) (1982) 447–455. doi:10.1016/0001-6160(82)90225-5 .3313] M. E. Gurtin, E. Fried, L. Anand, The Mechanics and Thermodynamics of Continua,Cambridge University Press, Cambridge, 2010.[14] H. Xiao, O. Bruhns, A. Meyers, Elastoplasticity beyond small deformations, ActaMechanica 182 (1-2) (2006) 31–111. doi:10.1007/s00707-005-0282-7 .[15] C. H. Rycroft, Y. Sui, E. Bouchbinder, An eulerian projection method for quasi-staticelastoplasticity, Journal of Computational Physics 300 (2015) 136–166. doi:10.1016/j.jcp.2015.06.046 .[16] W. Prager, An elementary discussion of definitions of stress rate, Quart. Appl. Math. 18(1960) 403–407.[17] C. Reina, S. Conti, Kinematic description of crystal plasticity in the finite kinematicframework: A micromechanical understanding of F = F e F p , Journal of the Mechanicsand Physics of Solids 67 (2014) 40–61. doi:10.1016/j.jmps.2014.01.014 .[18] C. Truesdell, Hypo-elasticity, Indiana Univ. Math. J. 4 (1955) 83–133.[19] R. Hill, A general theory of uniqueness and stability in elastic–plastic solids, Journal ofthe Mechanics and Physics of Solids 6 (3) (1958) 236–249. doi:10.1016/0022-5096(58)90029-2 .[20] J. B. Bell, P. Colella, H. M. Glaz, A second-order projection method for the incompressibleNavier–Stokes equations, Journal of Computational Physics 85 (2) (1989) 257–283. doi:10.1016/0021-9991(89)90151-4 .[21] A. Almgren, J. Bell, W. Szymczak, A numerical method for the incompressible Navier–Stokes equations based on an approximate projection, SIAM Journal on ScientificComputing 17 (2) (1996) 358–369. doi:10.1137/S1064827593244213 .[22] M. Sussman, A. S. Almgren, J. B. Bell, P. Colella, L. H. Howell, M. L. Welcome, Anadaptive level set approach for incompressible two-phase flows, J. Comput. Phys. 148 (1)(1999) 81–124. doi:10.1006/jcph.1998.6106 .[23] W. L. Briggs, V. E. Henson, S. F. McCormick, A Multigrid Tutorial: Second Edition,Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2000.[24] J. W. Demmel, Applied Numerical Linear Algebra, SIAM, 1997.[25] R. Courant, K. Friedrichs, H. Lewy, On the partial difference equations of mathematicalphysics, IBM Journal of Research and Development 11 (2) (1967) 215–234. doi:10.1147/rd.112.0215 .[26] B. Yang, M. L. Morrison, P. K. Liaw, R. A. Buchanan, G. Wang, C. T. Liu, M. Denda,Dynamic evolution of nanoscale shear bands in a bulk-metallic glass, Applied PhysicsLetters 86 (14) (2005) 141904. doi:10.1063/1.1891302 .3427] A. J. Chorin, A numerical method for solving incompressible viscous flow problems,Journal of Computational Physics 2 (1) (1967) 12–26. doi:10.1016/0021-9991(67)90037-X .[28] A. J. Chorin, Numerical solution of the Navier–Stokes equations, Mathematics ofComputation 22 (104) (1968) 745–762. doi:10.1090/S0025-5718-1968-0242392-2 .[29] M. L. Falk, J. S. Langer, Dynamics of viscoplastic deformation in amorphous solids,Phys. Rev. E 57 (6) (1998) 7192–7205. doi:10.1103/PhysRevE.57.7192 .[30] E. Bouchbinder, J. S. Langer, I. Procaccia, Athermal shear-transformation-zone theoryof amorphous plastic deformation. I. Basic principles, Phys. Rev. E 75 (3) (2007) 036107. doi:10.1103/PhysRevE.75.036107 .[31] J. S. Langer, Shear-transformation-zone theory of plastic deformation near the glasstransition, Phys. Rev. E 77 (2) (2008) 021502. doi:10.1103/PhysRevE.77.021502 .[32] E. Bouchbinder, J. S. Langer, Nonequilibrium thermodynamics of driven amorphousmaterials. I. Internal degrees of freedom and volume deformation, Phys. Rev. E 80 (3)(2009) 031131. doi:10.1103/PhysRevE.80.031131 .[33] J. J. Kruzic, Bulk metallic glasses as structural materials: A review, Advanced Engi-neering Materials 18 (8) (2016) 1308–1331. doi:10.1002/adem.201600066 .[34] T. C. Hufnagel, T. Jiao, Y. Li, L.-Q. Xing, K. T. Ramesh, Deformation and failureof Zr Ti Cu Ni Al bulk metallic glass under quasi-static and dynamic compression,Journal of Materials Research 17 (6) (2002) 1441–1445. doi:10.1557/JMR.2002.0214 .[35] M. L. Manning, J. S. Langer, J. M. Carlson, Strain localization in a shear transformationzone model for amorphous solids, Phys. Rev. E 76 (5) (2007) 056106. doi:10.1103/PhysRevE.76.056106 .[36] M. L. Manning, E. G. Daub, J. S. Langer, J. M. Carlson, Rate-dependent shear bandsin a shear-transformation-zone model of amorphous solids, Phys. Rev. E 79 (1) (2009)016110. doi:10.1103/PhysRevE.79.016110 .[37] B. A. Sun, S. Pauly, J. Hu, W. H. Wang, U. K¨uhn, J. Eckert, Origin of intermittentplastic flow and instability of shear band sliding in bulk metallic glasses, Phys. Rev. Lett.110 (2013) 225501. doi:10.1103/PhysRevLett.110.225501 .[38] J. Antonaglia, W. J. Wright, X. Gu, R. R. Byer, T. C. Hufnagel, M. LeBlanc, J. T. Uhl,K. A. Dahmen, Bulk metallic glasses deform via slip avalanches, Phys. Rev. Lett. 112(2014) 155501. doi:10.1103/PhysRevLett.112.155501 .[39] M. L. Falk, J. Langer, Deformation and failure of amorphous, solidlike materials,Annual Review of Condensed Matter Physics 2 (1) (2011) 353–373. doi:{10.1146/annurev-conmatphys-062910-140452} .3540] L. O. Eastgate, J. S. Langer, L. Pechenik, Dynamics of large-scale plastic deformationand the necking instability in amorphous solids, Phys. Rev. Lett. 90 (4) (2003) 045506. doi:10.1103/PhysRevLett.90.045506 .[41] C. H. Rycroft, F. Gibou, Simulations of a stretching bar using a plasticity model fromthe shear transformation zone theory, Journal of Computational Physics 231 (5) (2012)2155–2179. doi:10.1016/j.jcp.2011.10.009 .[42] A. Moriel, E. Bouchbinder, Necking instabilities in elastoviscoplastic materials, Phys.Rev. Materials 2 (2018) 073602. doi:10.1103/PhysRevMaterials.2.073602 .[43] C. H. Rycroft, E. Bouchbinder, Fracture toughness of metallic glasses: Annealing-inducedembrittlement, Phys. Rev. Lett. 109 (2012) 194301. doi:10.1103/PhysRevLett.109.194301 .[44] M. Vasoya, C. H. Rycroft, E. Bouchbinder, Notch fracture toughness of glasses: De-pendence on rate, age, and geometry, Phys. Rev. Applied 6 (2016) 024008. doi:10.1103/PhysRevApplied.6.024008 .[45] J. Ketkaew, W. Chen, H. Wang, A. Datye, M. Fan, G. Pereira, U. D. Schwarz, Z. Liu,R. Yamada, W. Dmowski, M. D. Shattuck, C. S. O’Hern, T. Egami, E. Bouchbinder,J. Schroers, Mechanical glass transition revealed by the fracture toughness of metallicglasses, Nature Communications 9 (1) (2018) 3271. doi:10.1038/s41467-018-05682-8 .URL https://doi.org/10.1038/s41467-018-05682-8 [46] A. S. Argon, Plastic deformation in metallic glasses, Acta Metallurgica 27 (1979) 47–58.[47] A. Wisitsorasak, P. G. Wolynes, Dynamical theory of shear bands in structural glasses,Proc. Natl. Acad. Sci. 114 (6) (2017) 1287–1292. doi:10.1073/pnas.1620399114 .[48] H. Xiao, O. T. Bruhns, A. Meyers, Hypo-elasticity model based upon the logarithmicstress rate, Journal of Elasticity 47 (1) (1997) 51–68. doi:10.1023/A:1007356925912 .[49] C. Truesdell, Hypoelastic shear, Journal of Applied Physics 27 (5) (1956) 441–447. doi:10.1063/1.1722399 .[50] J. Bardet, Finite element analysis of surface instability in hypo-elastic solids, ComputerMethods in Applied Mechanics and Engineering 78 (3) (1990) 273–296. doi:10.1016/0045-7825(90)90002-4 .[51] S. Hajarolasvadi, A. E. Elbanna, A new hybrid numerical scheme for modelling elasto-dynamics in unbounded media with near-source heterogeneities, Geophysical JournalInternational 211 (2) (2017) 851–864. doi:10.1093/gji/ggx337 .[52] X. Ma, A. Elbanna, Strain localization in dry sheared granular materials: A compactivity-based approach, Phys. Rev. E 98 (2018) 022906. doi:10.1103/PhysRevE.98.022906 .3653] T. J. Hughes, Numerical Implementation of Constitutive Models: Rate-IndependentDeviatoric Plasticity, Springer, 1984.[54] L. Anand, M. Kothari, A computational procedure for rate-independent crystal plasticity,Journal of the Mechanics and Physics of Solids 44 (4) (1996) 525–558. doi:10.1016/0022-5096(96)00001-4 .[55] G. Puglisi, L. Truskinovsky, Thermodynamics of rate-independent plasticity, Journalof the Mechanics and Physics of Solids 53 (3) (2005) 655–679. doi:10.1016/j.jmps.2004.08.004 .[56] J. Simo, R. Taylor, Consistent tangent operators for rate-independent elastoplasticity,Computer Methods in Applied Mechanics and Engineering 48 (1) (1985) 101–118. doi:10.1016/0045-7825(85)90070-2 .[57] R. Saye, Implicit mesh discontinuous galerkin methods and interfacial gauge methodsfor high-order accurate interface dynamics, with applications to surface tension dy-namics, rigid body fluidstructure interaction, and free surface flow: Part i, Journal ofComputational Physics 344 (2017) 647–682. doi:10.1016/j.jcp.2017.04.076 .[58] R. Saye, Implicit mesh discontinuous galerkin methods and interfacial gauge methodsfor high-order accurate interface dynamics, with applications to surface tension dy-namics, rigid body fluidstructure interaction, and free surface flow: Part ii, Journal ofComputational Physics 344 (2017) 683–723. doi:10.1016/j.jcp.2017.05.003 .[59] R. Saye, Interfacial gauge methods for incompressible fluid dynamics, Science Advances2 (6) (2016). doi:10.1126/sciadv.1501869 .[60] D. L. Brown, R. Cortez, M. L. Minion, Accurate projection methods for the incompressibleNavier–Stokes equations, J. Comput. Phys. 168 (2) (2001) 464–499. doi:10.1006/jcph.2001.6715 .[61] C. Min, F. Gibou, A second order accurate projection method for the incompressibleNavier–Stokes equations on non-graded adaptive grids, Journal of Computational Physics219 (2) (2006) 912–929. doi:10.1016/j.jcp.2006.07.019 .[62] Y. Zhang, A. L. Greer, Thickness of shear bands in metallic glasses, Applied PhysicsLetters 89 (7) (2006) 071907. doi:10.1063/1.2336598 .[63] A. Greer, Y. Cheng, E. Ma, Shear bands in metallic glasses, Materials Science andEngineering: R: Reports 74 (4) (2013) 71–132. doi:10.1016/j.mser.2013.04.001 .[64] C. A. Schuh, T. C. Hufnagel, U. Ramamurty, Mechanical behavior of amorphous alloys,Acta Materialia 55 (12) (2007) 4067–4109. doi:10.1016/j.actamat.2007.01.052 .[65] R. Maaß, J. F. L¨offler, Shear-band dynamics in metallic glasses, Advanced FunctionalMaterials 25 (16) (2015) 2353–2368. doi:10.1002/adfm.201404223 .3766] R. Maa, K. Samwer, W. Arnold, C. A. Volkert, A single shear band in a metallicglass: Local core and wide soft zone, Applied Physics Letters 105 (17) (2014) 171902. doi:10.1063/1.4900791 .[67] J. Lubliner, Plasticity Theory, Dover, New York, 2008.[68] D. L. Brown, Accuracy of projection methods for the incompressible navier-stokesequations (01 2003). doi:10.1142/9789812796837_0006 .[69] D. M. Summers, A. J. Chorin, Numerical vorticity creation based on impulse conservation,Proceedings of the National Academy of Sciences 93 (5) (1996) 1881–1885. doi:10.1073/pnas.93.5.1881 .[70] R. Cortez, An impulse-based approximation of fluid motion due to boundary forces,Journal of Computational Physics 123 (2) (1996) 341 – 353. doi:https://doi.org/10.1006/jcph.1996.0028 .[71] T. F. Buttke, Velicity Methods: Lagrangian Numerical Methods which Preserve theHamiltonian Structure of Incompressible Fluid Flow, 1993, pp. 39–57. doi:10.1007/978-94-015-8137-0_3 .[72] K. Kamrin, E. Bouchbinder, Two-temperature continuum thermomechanics of deformingamorphous solids, Journal of the Mechanics and Physics of Solids 73 (2014) 269–288. doi:10.1016/j.jmps.2014.09.009 .[73] E. Bouchbinder, J. S. Langer, Shear-transformation-zone theory of linear glassy dynamics,Phys. Rev. E 83 (2011) 061503. doi:10.1103/PhysRevE.83.061503 .[74] E. Bouchbinder, Effective temperature dynamics in an athermal amorphous plasticitytheory, Phys. Rev. E 77 (2008) 051505.[75] D. Loi, S. Mossa, L. F. Cugliandolo, Effective temperature of active matter, Phys. Rev.E 77 (2008) 051111. doi:10.1103/PhysRevE.77.051111 .[76] L. F. Cugliandolo, The effective temperature, Journal of Physics A: Mathematical andTheoretical 44 (48) (2011) 483001.[77] E. Bouchbinder, J. S. Langer, I. Procaccia, Athermal shear-transformation-zone theoryof amorphous plastic deformation. II.Analysis of simulated amorphous silicon, Phys. Rev.E 75 (3) (2007) 036108. doi:10.1103/PhysRevE.75.036108 .[78] M. T. Heath, Scientific Computing: An Introductory Survey, McGraw-Hill, 2002.[79] C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comp. Phys. 77 (2) (1988) 439–471. doi:10.1016/0021-9991(88)90177-5 . 3880] E. Gabriel, G. E. Fagg, G. Bosilca, T. Angskun, J. J. Dongarra, J. M. Squyres, V. Sahay,P. Kambadur, B. Barrett, A. Lumsdaine, R. H. Castain, D. J. Daniel, R. L. Graham, T. S.Woodall, Open MPI: Goals, concept, and design of a next generation MPI implementation,in: Proceedings, 11th European PVM/MPI Users’ Group Meeting, Budapest, Hungary,2004, pp. 97–104.[81] J. Xu, Iterative methods by space decomposition and subspace correction, SIAM Review34 (4) (1992) 581–613. doi:10.1137/1034116 .[82] M. Theillard, C. H. Rycroft, F. Gibou, A multigrid method on non-graded adaptiveoctree and quadtree cartesian grids, Journal of Scientific Computing 55 (1) (2013) 1–15. doi:10.1007/s10915-012-9619-2 .[83] J. Li, Z. L. Wang, T. C. Hufnagel, Characterization of nanometer-scale defects in metallicglasses by quantitative high-resolution transmission electron microscopy, Phys. Rev. B65 (2002) 144201. doi:10.1103/PhysRevB.65.144201 .[84] A. R. Hinkle, C. H. Rycroft, M. D. Shields, M. L. Falk, Coarse graining atomisticsimulations of plastically deforming amorphous solids, Phys. Rev. E 95 (2017) 053001. doi:10.1103/PhysRevE.95.053001 .[85] B. Sun, S. Pauly, J. Tan, M. Stoica, W. Wang, U. Khn, J. Eckert, Serrated flow andstickslip deformation dynamics in the presence of shear-band interactions for a zr-basedmetallic glass, Acta Materialia 60 (10) (2012) 4160–4171. doi:10.1016/j.actamat.2012.04.013doi:10.1016/j.actamat.2012.04.013