Parallel three-dimensional simulations of quasi-static elastoplastic solids

TL;DR

This paper introduces a parallel 3D simulation method for quasi-static hypo-elastoplastic solids, extending a projection approach based on incompressible Navier-Stokes equations, enabling detailed modeling of shear band formation in complex materials.

Contribution

It develops a novel 3D parallel geometric multigrid solver for quasi-static hypo-elastoplastic simulations, applicable to realistic materials and boundary conditions.

Findings

Successfully simulated shear band nucleation and growth in 3D.

Demonstrated applicability to bulk metallic glasses and heterogeneous materials.

Provided insights into shear band formation under various deformation conditions.

Abstract

Hypo-elastoplasticity is a flexible framework for modeling the mechanics of many hard materials under small elastic deformation and large plastic deformation. Under typical loading rates, most laboratory tests of these materials happen in the quasi-static limit, but there are few existing numerical methods tailor-made for this physical regime. In this work, we extend to three dimensions a recent projection method for simulating quasi-static hypo-elastoplastic materials. The method is based on a mathematical correspondence to the incompressible Navier-Stokes equations, where the projection method of Chorin (1968) is an established numerical technique. We develop and utilize a three-dimensional parallel geometric multigrid solver employed to solve a linear system for the quasi-static projection. Our method is tested through simulation of three-dimensional shear band nucleation and growth, a precursor to failure in many materials. As an example system, we employ a physical model of a bulk metallic glass based on the shear transformation zone theory, but the method can be applied to any elastoplasticity model. We consider several examples of three-dimensional shear banding, and examine shear band formation in physically realistic materials with heterogeneous initial conditions under both simple shear deformation and boundary conditions inspired by friction welding.

Figures

Tables

Table 1. Table 1 : Material parameters used in this study, for both linear elasticity and the STZ model of amorphous plasticity. The Boltzmann constant k B subscript 𝑘 𝐵 k_{B} is used to convert energetic values to temperatures.

Parameter	Value
Young’s modulus $E$	101 GPa
Poisson ratio $\nu$	0.35
Bulk modulus $K$	122 GPa
Shear modulus $\mu$	37.4 GPa
Density ${\rho }_0$	$6125{\text{ kg m}}^{-3}$
Shear wave speed $c_s$	$2.47{\text{ km s}}^{-1}$
Yield stress $s_Y$	0.85 GPa
Molecular vibration timescale ${\tau }_0$	${10}^{-13}$ s
Typical local strain ${ϵ}_0$	0.3
Effective heat capacity $c_0$	0.4
Typical activation barrier $\Delta /k_B$	8000 K
Typical activation volume $\Omega$	300 Å3
Thermodynamic bath temperature $T$	400 K
Steady state effective temperature ${\chi }_{\infty }$	900 K
STZ formation energy $e_z/k_B$	21000 K

Table 2. Table 2 : Data describing the total time per simulation, the time spent in all V-cycles per simulation, the total number of V-cycles, the average time per V-cycle, and the weak scaling ratios for total time and V-cycle time. Values of N 𝑁 N and p 𝑝 p were chosen to keep the number of grid points per processor N 3 / ( 2 ​ p ) superscript 𝑁 3 2 𝑝 N^{3}/(2p) approximately constant. The weak scaling (WS) ratio is computed using the formula w = r ​ t 1 / t p 𝑤 𝑟 subscript 𝑡 1 subscript 𝑡 𝑝 w=rt_{1}/t_{p} where t 1 subscript 𝑡 1 t_{1} is the time for the simulation with one processor, t p subscript 𝑡 𝑝 t_{p} is the time for the simulation with p 𝑝 p processors, and r = N p 3 / ( p ​ N 1 3 ) 𝑟 subscript superscript 𝑁 3 𝑝 𝑝 superscript subscript 𝑁 1 3 r=N^{3}_{p}/(pN_{1}^{3}) is a small adjustment factor for deviations from perfectly equal work.

	$N=256$ $p=32$	$N=204$ $p=16$	$N=164$ $p=8$	$N=128$ $p=4$	$N=102$ $p=2$	$N=80$ $p=1$
Total time (hours)	3.69	2.31	2.36	1.94	1.74	1.52
V-cycle time (hours)	2.75	1.71	1.74	1.37	1.20	1.00
# of V-cycles	17939	17539	17171	16511	15575	15073
Time/V-cycle (seconds)	0.55	0.35	0.37	0.30	0.28	0.24
WS ratio (total time)	0.42	0.68	0.69	0.80	0.90	1.0
WS ratio (V-cycle time)	0.37	0.61	0.62	0.75	0.86	1.0

Table 3. Table 3 : Data describing the total time per simulation, the time spent in all V-cycles per simulation, the total number of V-cycles, and the average time per V-cycle. p 𝑝 p is fixed at 32 while N 𝑁 N is varied.

	$N=128$	$N=256$	$N=384$	$N=512$	$N=768$
Total time (hours)	1.61	6.12	22.44	42.86	105.97
V-cycle time (hours)	1.27	3.89	14.28	30.23	75.90
# of V-cycles	20480	18367	24675	28175	19689
Time/V-cycle (seconds)	0.22	0.76	2.08	3.87	13.88

Equations

$$\frac{\mathcal{D}\mathbf{\sigma }}{\mathcal{D}t}=\mathbf{C}:{\mathbf{D}}^{\text{el}}=\mathbf{C}:\left(\mathbf{D}-{\mathbf{D}}^{\text{pl}}\right).$$

$$\rho \frac{d\mathbf{v}}{dt}=\nabla \cdot \mathbf{\sigma }$$

$$\nabla \cdot \mathbf{\sigma }\approx 0.$$

$$\frac{d\mathbf{v}}{dt}=-\nabla p+\nu {\nabla }^2\mathbf{v}.$$

$$\frac{d\rho }{dt}=-\rho \left(\nabla \cdot \mathbf{v}\right),$$

$$\nabla \cdot \mathbf{v}=0.$$

$${\mathbf{v}}^{\ast }=\mathbf{v}+\nabla \varphi ,$$

$${\mathbf{v}}_t^{\ast }+\left(\mathbf{v}\cdot \nabla \right)\mathbf{v}+\nabla q=\nu {\nabla }^2{\mathbf{v}}^{\ast },$$

$$\nabla p=\nabla \left(q+{\varphi }_t\right)-\nu {\nabla }^2\nabla \varphi ,$$

$${\nabla }^2\varphi =\nabla \cdot {\mathbf{v}}^{\ast }.$$

$$\mathbf{v}={\mathbf{v}}^{\ast }-\nabla \varphi .$$

$$⟨\mathbf{v},\mathbf{u}⟩={\int }_{\Omega }\mathbf{v}\cdot \mathbf{u}{d}^3\mathbf{x},$$

$$⟨{\mathbf{v}}^{n+1}-{\mathbf{v}}^n,{\mathbf{v}}^{n+1}-{\mathbf{v}}^{\ast }⟩$$

$$={\int }_{\Omega }\left(\nabla \cdot {\mathbf{v}}^{n+1}-\nabla \cdot {\mathbf{v}}^n\right)\varphi (\mathbf{x})d^3\mathbf{x}=0,$$

$${\mathbf{\sigma }}^{\ast }=\mathbf{\sigma }+\mathbf{C}:\nabla \mathbf{\Phi },$$

$${\mathbf{\sigma }}_t^{\ast }+\left(\mathbf{v}\cdot \nabla \right)\mathbf{\sigma }=\mathbf{L}\mathbf{\sigma }+\mathbf{\sigma }{\mathbf{L}}^{\mathsf{T}}-\text{Tr}(\mathbf{L})\mathbf{\sigma }+\mathbf{C}:\left(\nabla \mathbf{q}-{\mathbf{D}}^{\text{pl}}\right).$$

$$\mathbf{C}:\mathbf{D}=\mathbf{C}:\nabla \left(\mathbf{q}-{\mathbf{\Phi }}_t\right),$$

$$\nabla \cdot \left(\mathbf{C}:\nabla \mathbf{\Phi }\right)=\nabla \cdot {\mathbf{\sigma }}^{\ast }.$$

$$⟨\mathbf{\sigma },{\mathbf{\sigma }}^{\prime }⟩={\int }_{\Omega }\mathbf{\sigma }:\mathbf{S}:{\mathbf{\sigma }}^{\prime }{d}^3\mathbf{x},$$

$$⟨{\mathbf{\sigma }}^{n+1}-{\mathbf{\sigma }}^n,{\mathbf{\sigma }}^{n+1}-{\mathbf{\sigma }}^{\ast }⟩$$

$$=\int \left({\mathbf{\sigma }}^{n+1}-{\mathbf{\sigma }}^n\right):\nabla \mathbf{\Phi }{d}^3\mathbf{x}$$

$$=-\int \left(\nabla \cdot {\mathbf{\sigma }}^{n+1}-\nabla \cdot {\mathbf{\sigma }}^n\right)\cdot \mathbf{\Phi }{d}^3\mathbf{x}=0.$$

$${\tau }_0{D}^{\text{pl}}=e^{-e_z/k_B\chi }\mathcal{C}(\bar{s},T)\left(1-\frac{s_Y}{\bar{s}}\right),$$

$$\mathcal{R}(\pm \bar{s},T)=\exp \left(-\frac{\Delta \mp \Omega {ϵ}_0\bar{s}}{k_{B}T}\right).$$

$$\mathcal{C}(\bar{s},T)=e^{-\Delta /k_{B}T}\cosh \left(\frac{\Omega {ϵ}_0\bar{s}}{k_{B}T}\right).$$

$$c_0\frac{d\chi }{dt}=\frac{\left({\mathbf{D}}^{\text{pl}}:{\mathbf{\sigma }}_0\right)}{s_Y}\left({\chi }_{\infty }-\chi \right)+l^2\nabla \cdot \left(D^{\text{pl}}\nabla \chi \right).$$

$$\rho \frac{{\mathbf{v}}^{n+1}-{\mathbf{v}}^n}{\Delta t}$$

$$\frac{{\mathbf{\sigma }}^{n+1}-{\mathbf{\sigma }}^n}{\Delta t}$$

$$+\text{Tr}({\mathbf{L}}^n){\mathbf{\sigma }}^n+\mathbf{C}:\left({\mathbf{D}}^n-\frac{{\left(D^{\text{pl}}\right)}^n}{{\bar{s}}^n}{\mathbf{\sigma }}_0^n\right),$$

$$c_0\frac{{\chi }^{n+1}-{\chi }^n}{\Delta t}$$

$${\mathbf{\sigma }}^{\ast }={\mathbf{\sigma }}^n+\Delta t\left({\mathbf{L}}^n{\mathbf{\sigma }}^n+{\mathbf{\sigma }}^n{\left({\mathbf{L}}^{\mathsf{T}}\right)}^n-\text{Tr}({\mathbf{L}}^n){\mathbf{\sigma }}^n-\mathbf{C}:\left(\frac{D^{\text{pl}}}{{\bar{s}}^n}{\mathbf{\sigma }}_0^n\right)\right).$$

$${\mathbf{D}}^{n+1}=\frac{1}{2}\left(\left(\nabla {\mathbf{v}}^{n+1}\right)+{\left(\nabla {\mathbf{v}}^{n+1}\right)}^{\mathsf{T}}\right)$$

$${\mathbf{\sigma }}^{n+1}={\mathbf{\sigma }}^{\ast }+\Delta t\left(\mathbf{C}:{\mathbf{D}}^{n+1}\right).$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Modeling the 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 localization of strain may be driven by geometrical effects such as necking [2, 3, 4, 5], or by instabilities in the material response, such as due to plastic yielding. In the latter case, positive feedback causes regions that have already plastically deformed to yield further, often resulting in a localized plane of strain called a shear band [6, 7, 8, 9, 10]. Shear bands have been studied analytically [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, and herein we develop a three-dimensional extension of a two-dimensional algorithm [13] suitable for simulating deformation in hard elastoplastic materials. As one application, we demonstrate the ability of the method to capture fine-scale features of shear banding in the deformation of amorphous materials.

For stress levels below the material yield stress, an elastoplastic material deforms purely elastically and returns to its undeformed state upon removal of the load. When the yield stress is reached, the material begins to deform plastically, leading to permanent, irreversible deformation that persists beyond load removal [14]. Elastoplasticity is complex, and it admits a number of mathematical descriptions [15, 13], each of which amounts to a specification of the interaction of the elastic and plastic components of deformation at a microscopic level [16, 17]. In this paper, we utilize the hypo-elastoplasticity model [18], in which the Eulerian rate of deformation tensor is decomposed additively into elastic and plastic parts, $\mathbf{D}=\mathbf{D}^{\text{el}}+\mathbf{D}^{\text{pl}}$ [19]. Hypo-elastoplasticity has some drawbacks: for example, it is generally only used for linear elastic deformation, because it is difficult to represent nonlinear elasticity through $\mathbf{D}^{\text{el}}$ alone [20, 21]. The formalism is also based on the material velocity rather than displacement field, and hence can lead to accumulation of numerical integration errors [22, 23]. Despite these flaws, hypo-elastoplasticity is well-suited to problems with small elastic deformation and large plastic deformation, which makes it particularly appropriate for studying failure and shear banding in hard materials.

On the other hand, the hypo-elastoplastic formulation has several numerical advantages. Since it is based on the Eulerian rate of deformation tensor, it is well-suited for a fixed-grid framework. Fixed grids have simpler topologies than their Lagrangian counterparts, and are easier to program and parallelize. This is particularly important in three dimensions, where the computational expense mandates parallelization, and where the implementation difficulty increases relative to lower dimensional simulations. Fixed-grid methods are also the methods of choice for fluid simulation [24, 25, 26], and they allow a wider range of numerical linear algebra techniques to be used, such as the geometric multigrid method [27, 28].

The additive decomposition of $\mathbf{D}$, coupled with the linear-elastic constitutive relation and a continuum formulation of Newton’s second law, leads to a closed hyperbolic system of partial differential equations (PDEs) for the material velocity, stress, and variables intrinsic to the plasticity model. To properly resolve elastic waves, the timestep $\Delta t$ for an explicit numerical scheme is restricted by the well-known Courant–Friedrichs–Lewy (CFL) condition [29]. The CFL condition states that $\Delta t$ must obey a constraint for numerical stability: $\Delta t\leq\frac{h}{c_{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 elastic wave travel time [30]. The CFL condition thus poses a prohibitive limit on the timestep for probing realistic timescales and system sizes, and the development of alternative simulation approaches which avoid resolving elastic waves is necessary. By looking at the limit of long times and small velocities, one can show that Newton’s second law can be replaced by a constraint that the stresses must remain in quasi-static equilibrium. However, these equations are no longer a hyperbolic system, since the ability to time-integrate the velocity field is lost.

Recently, Rycroft et al. [13] demonstrated a mathematical analogy between quasi-static elastoplasticity and the incompressible Navier–Stokes equations for fluid flow. The incompressible Navier–Stokes equations combine an equation for the velocity – dependent on the fluid pressure – with a requirement that the velocity must be divergence free. Much like in the case of quasi-static elastoplasticity, the divergence-free requirement on the velocity field is obtained as a limit of a PDE for the pressure. A well-known algorithm in this setting is the projection method of Chorin [31, 32]. In Chorin’s method, an intermediate step is taken for the velocity field, but this intermediate velocity does not obey the incompressibility constraint. An elliptic problem is solved that simultaneously enforces the incompressibility constraint and enables computation of the pressure. Rycroft et al. translated Chorin’s projection method to quasi-static elastoplasticity by exploiting the mathematical analogy between the two sets of limiting equations. The new method was quantitatively tested on some simple two-dimensional systems using a multi-threaded single-processor implementation. It was shown to be practical and efficient, enabling the use of timesteps that are orders of magnitude larger than those required by the CFL condition.

Analogies between elastoplastic deformation and incompressible flow have been noted previously, and have been used to employ Vanka-type smoothers for solid mechanics problems [33, 34]. However, there are several important differences between this prior work and our own: we employ a formalism based on the material velocity field, while previous approaches solve for the displacement field. Similarly, to our knowledge, our work is the first projection-type method for elastoplasticity based on Chorin’s method.

In a similar vein, our method has strong similarity to numerical algorithms for low-Mach flows. Such methods avoid the resolution of acoustic waves yet retain compressible effects in the fluid, and much like incompressible fluid dynamics, result in an elliptic equation for the pressure. They have been used to study Type Ia supernovae [35] and nuclear combustion [36], and have their roots in an earlier approach due to Majda and Sethian [37]. Analytical justification for the limiting procedures taken has been provided by Klainerman and Majda [38, 39]. It is possible that similar approaches could be developed for alternative limits of the equations of hypo-elastoplasticity.

In this paper, we extend the quasi-static projection method to three dimensions. Due to the greatly increased problem size, we develop a distributed parallel implementation, which creates algorithmic challenges. In particular, the projection method requires that we solve a coupled elliptic equation for the components of velocity. We develop a custom parallel geometric multigrid code to solve the resulting linear system. As a physical testbed for our methodology, we employ an athermal formulation of the shear transformation zone (STZ) theory of Falk, Langer, Bouchbinder and coworkers as a plasticity model for a bulk metallic glass [40, 41, 42, 43]. Metallic glasses naturally lend themselves to study through the hypo-elastoplasticity framework, as their elastic deformation is generally small and well-described by a linear theory, yet they can exhibit significant plastic deformation [44]. Their elastic moduli are typically on the order of 10–100 GPa, and hence experimental loading conditions often place samples in the quasi-static regime [45]. They also present interesting and poorly understood fundamental physics [46, 47, 48, 49], which can be difficult to probe in dimensions greater than one without the numerical methods presented here.

Simulation studies of bulk metallic glasses have been essential to our understanding of their properties. The development of the original STZ theory was guided by observations of molecular dynamics simulations [50]. Simulations of the necking instability in a bar under uniaxial tension using the STZ plasticity model highlight the interplay between elastic and plastic deformation [51, 52, 53]. The physical mechanisms of fracture in BMGs were explored using the two-dimensional projection method [54], subsequently allowing the fracture toughness of BMGs to be predicted across a wide range of experimental conditions [55]. Later experimental measurements due to Ketkaew et al. demonstrated that these simulation-based predictions were quantitatively correct [56]. Indeed, testable predictions for complex amorphous systems such as BMGs are rare, and the development of efficient numerical methods such as the ones presented here provide a way to generate them and to guide future experimental inquiry.

The STZ theory is a useful test case for our method both physically and numerically, but the numerical methodology is general and can be used for many plasticity models within the hypo-elastoplasticity framework. These could include free-volume based models of BMGs [57], plasticity models based on the random first-order transition theory of the glass transition [58], hypo-elastic materials [59, 60, 61], geophysical models [62, 63], and rate-independent plasticity models [64, 65, 66, 67]. We also emphasize that Chorin’s projection method represents a first step towards more sophisticated projection-based algorithms such as gauge methods [68, 69, 70] and pressure-Poisson methods [71, 72], and that we have laid the groundwork here to generalize these algorithms to the case of hypo-elastoplasticity.

Under loading, BMGs exhibit shear bands [73, 74], which rapidly lead to material failure [30, 75, 76] and are one of the primary limitations in employing BMGs in applications [77]. Analytical work probing shear bands in amorphous materials is difficult, particularly in two or three dimensions, which highlights a need for computational investigations. The development of our method enables the study of shear-banding in three dimensions at large scale and high resolution without excessive computational expense. Our simulations demonstrate that this scale and resolution is indeed necessary, and expose interesting fine-scale and uniquely three-dimensional features of shear banding. Our methodology opens the door to future studies probing the shape, structure, and topology of shear bands, as well as the mechanism and statistical properties of their formation.

The structure of this paper is as follows. In Section 2, we describe the relation between Chorin’s projection method and the projection algorithm for hypo-elastoplasticity employed here. 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 equations in the non-quasi-static limit. We also discuss our parallel multigrid implementation used for the stress projection. In Section 4, we demonstrate convergence between the explicit and projection methods in a regime in which the two are expected to produce similar results, and study several interesting examples of shear banding dynamics in a metallic glass.

2 Projection methods for fluid dynamics and hypo-elastoplasticity

2.1 Hypo-elastoplasticity

We consider an elastoplastic material with Cauchy stress tensor $\boldsymbol{\sigma}(\mathbf{x},t)$ and velocity field $\mathbf{v}(\mathbf{x},t)$ at a position $\mathbf{x}$ and time $t$. The total rate of deformation tensor $\mathbf{D}$ is defined as the symmetric part of the velocity gradient, $\mathbf{D}=\frac{1}{2}(\mathbf{L}+\mathbf{L}^{\mathsf{T}})$ with $\mathbf{L}=\nabla\mathbf{v}$. For any field $f(\mathbf{x},t)$, we define the advective time derivative by $\frac{df}{dt}=\frac{\partial f}{\partial t}+\left(\mathbf{v}\cdot\nabla\right){f}$. The fundamental assumption of hypo-elastoplasticity is that the rate of deformation tensor can be additively decomposed into a sum of elastic and plastic parts, $\mathbf{D}=\mathbf{D}^{\text{el}}+\mathbf{D}^{\text{pl}}$.

For stiff elastoplastic materials with small elastic deformation, the linear elastic constitutive law provides an accurate description,

[TABLE]

In Eq. 1, $\mathbf{C}$ is the fourth-rank stiffness tensor, taken to be homogeneous and isotropic. With Lamé’s first parameter $\lambda$ and shear modulus $\mu$, the components of $\mathbf{C}$ are given by $C_{ijkl}=\lambda\delta_{ij}\delta_{kl}+\mu\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right)$ [78]. The time derivative $\frac{\mathcal{D}\boldsymbol{\sigma}}{\mathcal{D}t}=\frac{d\boldsymbol{\sigma}}{dt}-\mathbf{L}\boldsymbol{\sigma}-\boldsymbol{\sigma}\mathbf{L}^{T}+\text{Tr}(\mathbf{L})\boldsymbol{\sigma}$ denotes the Truesdell objective stress rate.

From Newton’s second law, the material velocity obeys the equation

[TABLE]

where $\rho$ denotes the material density. Equations 1 & 2 form a hyperbolic system of equations for the stress and velocity fields, which can be solved explicitly using standard finite-difference simulation methods. This hyperbolic system will resolve elastic waves, and so the timestep $\Delta t$ and grid spacing $\Delta x$ must satisfy the CFL condition $\Delta t\leq\Delta x/c_{e}$ for numerical stability, where $c_{e}$ is an elastic wave speed. In materials such as metals and metallic glasses, elastic waves travel on the order of kilometers per second. Spatial discretizations capable of resolving fine-scale features of interesting physical phenomena in these materials can be as small as micrometers. For $\Delta x=1\text{~{}\textmu{}m}$ and $c_{e}=1\text{~{}km/s}$, the CFL condition requires $\Delta t\leq 1\text{~{}ns}$, an extreme restriction for phenomena that occur on realistic timescales of hours or days.

2.2 Quasi-static hypo-elastoplasticity

We consider a scenario in which plastic deformation occurs on a timescale much greater than the time for waves to propagate through the material. In this setting, macroscopic plastic deformation takes place due to the accumulation of small velocity gradients over long times. The details of a limiting procedure describing this physical regime were performed in previous two-dimensional work [13] and will not be reproduced here.

In this quasi-static limit, the equation for the velocity in Eq. 2 can be approximately replaced by a constraint on the stress

[TABLE]

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 evolution equation for the stress in Eq. 1 is unaffected by the limiting procedure, and hence Eq. 1 must be 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 $\mathbf{v}$ appears in Eq. 1 through $\mathbf{D}$, but there is no longer an equation that can be integrated explicitly to solve for it. It is also not guaranteed that solutions of Eq. 1 subject to the constraint in Eq. 3 will agree with solutions of Eq. 1 and Eq. 2.

2.3 Incompressible fluid dynamics

We now demonstrate an analogy between the computational issues presented in the previous section and those encountered in incompressible fluid dynamics. Consider a fluid with velocity $\mathbf{v}$, pressure $p$, and density $\rho$. The fluid velocity field obeys the Navier–Stokes equation,

[TABLE]

If the fluid is compressible, then its density satisfies

[TABLE]

along with an equation of state linking the fluid density to the fluid pressure. Using an explicit scheme to solve the hyperbolic system in Eqs. 4 & 5 will resolve sound waves in the fluid, 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,

[TABLE]

This limit reduces the coupled partial differential equations for the pressure and velocity to a single constrained equation for the velocity. The pressure is present in the equation for the velocity, though its evolution equation has been exchanged for the incompressibility constraint; this is much like the quasi-static limit of hypo-elastoplasticity described in the preceding section.

Chorin [31, 32] developed a numerical method for this system of equations that involves the use of an orthogonal projection, spurring significant research into related algorithms [71, 72]. Such projection methods proceed via a two-step procedure, where an intermediate velocity $\mathbf{v}^{*}$ is first computed which does not obey the incompressibility constraint. $\mathbf{v}^{*}$ is then orthogonally projected onto the manifold of divergence-free solutions through the solution of an elliptic problem for an auxiliary field related to the pressure. The process of projection simultaneously enforces the constraint and enables computation of the pressure field.

One typical approach is to employ a Hodge decomposition [71, 72],

[TABLE]

where $\mathbf{v}$ is the desired divergence-free velocity field and $\phi$ is an auxiliary field. One then updates $\mathbf{v}^{*}$ for a fixed interval of time via the equation

[TABLE]

where $\nabla q$ is an approximation to the pressure gradient. Substituting Eq. 7 into Eq. 8 leads to a formula for the pressure

[TABLE]

from which $p$ can be computed. The divergence of Eq. 7 implies that $\nabla\cdot\mathbf{v}=0$ if $\phi$ is such that

[TABLE]

The Poisson problem in Eq. 10 can be solved for $\phi$ using standard techniques of numerical linear algebra, and the projection can be completed by computing

[TABLE]

Boundary conditions on $\phi$ in Eq. 10 depend on the physical scenario of interest, and are critical for obtaining higher-order methods [71, 72]. The algorithm typically proceeds by advancing Eq. 8 for a single timestep, computing $\mathbf{v}$ according to Eq. 11, and then setting $\mathbf{v}^{*}=\mathbf{v}$ for the beginning of the next timestep [79]. This procedure is schematically represented in discrete-time in Fig. 1(a). Projection methods avoid the CFL condition associated with compressive waves in the fluid, and hence can use significantly larger timesteps than explicit methods. An extension of projection methods known as gauge methods do not reset $\mathbf{v}^{*}=\mathbf{v}$ at the end of timestep, and instead allow it to continue to evolve during the computation [68, 69, 80, 81, 82]

It is possible to demonstrate that solving Eq. 10 represents an orthogonal projection. We define the inner product between two vector-valued fields,

[TABLE]

where $\Omega$ is the simulation domain. Using this inner product, we can compute

[TABLE]

thereby establishing that the projection $\mathbf{v}^{n+1}-\mathbf{v}^{*}$ is orthogonal to the difference between the two velocity fields, $\mathbf{v}^{n+1}-\mathbf{v}^{n}$.

2.4 A family of projection methods for hypo-elastoplasticity

We now formulate a three-dimensional projection method for solving Eq. 1 subject to the quasi-static constraint Eq. 3. We define an intermediate stress

[TABLE]

where $\boldsymbol{\Phi}(\mathbf{x},t)$ is an auxiliary vector field. We can solve for $\boldsymbol{\sigma}^{*}$ by dropping the $\mathbf{C}:\mathbf{D}$ term in Eq. 1,

[TABLE]

In Eq. 15, $\mathbf{q}$ represents an approximation to the velocity $\mathbf{v}$. Substituting Eq. 14 into Eq. 15, we find

[TABLE]

from which $\mathbf{D}$ can be computed. Taking the divergence of Eq. 14 and requiring $\nabla\cdot\boldsymbol{\sigma}=\mathbf{0}$, $\boldsymbol{\Phi}$ must satisfy the equation

[TABLE]

Equation 17 is a linear system with source term $\nabla\cdot\boldsymbol{\sigma}^{*}$ that can be solved for $\boldsymbol{\Phi}$. Hence, we can devise a projection method for quasi-static hypo-elastoplasticity by evolving Eq. 15 for a single timestep, finding $\boldsymbol{\Phi}$ according to Eq. 17, and projecting $\boldsymbol{\sigma}^{*}$ onto the manifold of divergence-free solutions by computing $\boldsymbol{\sigma}=\boldsymbol{\sigma}^{*}-\mathbf{C}:\nabla\boldsymbol{\Phi}$. The algorithm may then proceed by setting $\boldsymbol{\sigma}^{*}=\boldsymbol{\sigma}$ for the start of the next timestep. We represent this algorithm schematically in Fig. 1(b). A projection method is defined by the choice of the approximate velocity field $\mathbf{q}$, the auxiliary vector field $\boldsymbol{\Phi}$, and the integration method for Eq. 15. By instead allowing $\boldsymbol{\sigma}^{*}$ to evolve over the course of the computation rather than setting $\boldsymbol{\sigma}^{*}=\boldsymbol{\sigma}$ after each timestep, we expect it to be possible to develop gauge methods for quasi-static hypo-elastoplasticity in an analogous manner to the case of fluid dynamics.

As in the case of fluid dynamics, we can show that the projection is orthogonal in a suitable inner product. To do so, we define an inner product between two stress tensors as in [13],

[TABLE]

where $\mathbf{S}=\mathbf{C}^{-1}$ is the stiffness tensor. Equation 18 computes the elastic strain energy of a material with stress field $\boldsymbol{\sigma}$ and strain field $\mathbf{S}:\boldsymbol{\sigma}^{\prime}$, or vice-versa by symmetry. Because $\mathbf{S}$ is a symmetric positive definite tensor for physically realistic Lamé parameters, this definition is an inner product. By explicit computation,

[TABLE]

3 Numerical implementation

In this section, we describe an implementation of an explicit forward Euler method to solve 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 STZ theory.

3.1 Plasticity model

As a plasticity model for a metallic glass, we use an athermal form of the STZ theory suitable for studying diverse materials including BMGs below the glass transition temperature, dense granular materials, and soft materials such as foams or colloidal glasses [41, 42]. Within the STZ theory, irreversible molecular rearrangements are assumed to occur sporadically throughout an otherwise elastic material, and each rearrangement induces a small increment of strain. The accumulation of many such events leads to macroscopic plastic deformation. It is assumed that when local stresses surpass the material yield stress $s_{y}$, these rearrangements occur at rare and localized sites known as STZs. Thermal fluctuations of the atomic configuration are neglected in the athermal formulation: molecular rearrangements are entirely driven by external mechanical forces. Thermal theories introduce additional coupling between a configurational subsystem governing the rearrangements that occur at the STZs, and a kinetic/vibrational subsystem governing the thermal vibrations of atoms in their cage of nearest neighbors [83].

STZs may be conceptualized as clusters of atoms predisposed to configurational rearrangements when subjected to external shear [41]. Each rearrangement corresponds to a transition in the configurational energy landscape; these transitions are usually towards a lower-energy configuration, but there is a small probability for a reverse transition. Before the application of external shear, the material sample is at a local minimum. External shear alters the shape of 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 disorder temperature denoted by $\chi$ [84, 85, 86, 87]. $\chi$ governs the out-of-equilibrium configurational degrees of freedom of the material and has many properties of the usual temperature: it is measured in Kelvin, and can be obtained as the derivative of a configurational energy with respect to a configurational entropy [50]. $\chi$ is distinct from the thermodynamic temperature $T$, though it plays the same role for the configurational subsystem as $T$ does for the kinetic/vibrational subsystem.

The plastic rate of deformation tensor is proportional to the deviatoric part of the stress tensor $\boldsymbol{\sigma}_{0}=\boldsymbol{\sigma}-\frac{1}{3}\mathbf{I}\operatorname{tr}(\boldsymbol{\sigma})$, so that $\mathbf{D}^{\text{pl}}=D^{\text{pl}}\frac{\boldsymbol{\sigma}_{0}}{\bar{s}}$. $\bar{s}$ is a local stress measure essentially given by the Frobenius norm of the deviatoric stress tensor, $\bar{s}^{2}=\frac{1}{2}\sum_{ij}\sigma_{0,ij}^{2}$. The magnitude of the plastic rate of deformation is given by

[TABLE]

where $\tau_{0}$ is a molecular vibration timescale, $e_{z}$ is a typical STZ formation energy, and $k_{B}$ is the Boltzmann constant. $\mathcal{C}(\bar{s},T)$ represents the total STZ transition rate. With $\mathcal{R}(\pm\bar{s},T)$ denoting the forward and reverse rates between two configurational states, the total transition rate is $\mathcal{C}(\bar{s},T)=\frac{1}{2}\left(\mathcal{R}(\bar{s},T)+\mathcal{R}(-\bar{s},T)\right)$. The transitions follow a linearly stress-biased thermal activation process,

[TABLE]

$\Delta$ is a typical energetic barrier for a transition, $\Omega$ is a typical STZ volume, and $\epsilon_{0}$ is a typical local strain due to an STZ transition. While thermal fluctuations are neglected in the athermal model, the thermodynamic temperature still sets the magnitude of transition rates in the system. Using the form Eq. 21 yields the overall transition rate

[TABLE]

The effective temperature satisfies [41, 88, 46, 47]

[TABLE]

Equation 23 consists of a term causing growth to an asymptotic value $\chi_{\infty}$ and a diffusive term with length scale $l$. Both saturation to $\chi_{\infty}$ and diffusion occur in response to plastic deformation. The term $\mathbf{D}^{\text{pl}}:\boldsymbol{\sigma}_{0}$ is the energy dissipation rate due to externally applied mechanical work – STZs are created and annihilated proportionally – and $c_{0}$ is an effective heat capacity. Eq. 23 is thus essentially a heat equation, representing the first law of thermodynamics for the configurational subsystem [41]. The interdependence of Eqs. 20 & 23 enables the development of shear bands via positive feedback, as increasing $\chi$ also increases $D^{\text{pl}}$ [46, 47].

3.2 Explicit method

For the explicit method, we discretize Eqs. 1, 2, and 23 with a forward Euler step. This leads to the coupled set of discrete-time equations

[TABLE]

The small viscous stress term $\kappa\nabla^{2}\mathbf{v}$ in Eq. 24 is artificially imposed for numerical stability, but it is not needed in the quasi-static method. This term is physically reasonable, and could represent damping within the material. The requirement of damping in the explicit scheme for stability should be expected. If one focuses on elastic waves in the linearized regime then Eq. 26 can be neglected, and the leading-order terms in Eqs. 24 & 25 are $\nabla\cdot\boldsymbol{\sigma}^{n}$ and $\mathbf{C}:\mathbf{D}^{n}$, respectively. Hence Eqs. 24 & 25 form a first-order hyperbolic system for stress and velocity [89, 52]. We implement the explicit method using centered differencing for the spatial derivatives (Sec. 3.4), and it is well-known that such discretizations are unstable in the absence of any damping [90]. The scheme can be stabilized with a damping $\kappa=\kappa^{\prime}h$ where $\kappa^{\prime}$ is a constant that we determine empirically.

In three dimensions, the damping term induces a restriction on the timestep $\Delta t\leq\frac{h^{2}}{6\kappa}$. Hence, if $\kappa$ is viewed as a physical constant, this condition is more restrictive than the CFL condition. Here, we choose $\kappa$ to scale with the grid spacing as required for stability, and hence the timestep restriction scales in the same way as the CFL restriction.

3.3 Quasi-static method

We now formulate a specific three-dimensional projection method for solving Eq. 1 subject to the quasi-static constraint Eq. 3. We first neglect the $\mathbf{C}:\mathbf{D}$ term in Eq. 1 and compute an intermediate stress $\boldsymbol{\sigma}^{*}$,

[TABLE]

If the velocity at the next timestep $\mathbf{v}^{n+1}$ were known, we could compute

[TABLE]

and complete the forward Euler step in Eq. 27 as

[TABLE]

Taking the divergence of Eq. 29, requiring $\nabla\cdot\boldsymbol{\sigma}^{n+1}=0$, and rearranging terms leads to the equation

[TABLE]

Equation 30 is a linear system for the velocity $\mathbf{v}^{n+1}$ involving mixed spatial derivatives. The source term is given by $-\nabla\cdot\boldsymbol{\sigma}^{*}$. After solving for $\mathbf{v}^{n+1}$, it can be used to compute $\boldsymbol{\sigma}^{n+1}$ via Eq. 29. Through this process, $\boldsymbol{\sigma}^{*}$ is projected down to the manifold of divergence-free solutions to arrive at $\boldsymbol{\sigma}^{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 can nevertheless be solved rapidly via standard techniques of numerical linear algebra such as the multigrid method. The multigrid method relies on the Gauss–Seidel method for iterative smoothing of the solution, and Gauss–Seidel smoothing is guaranteed to converge if either the linear system is (A) weakly diagonally dominant, or (B) symmetric positive definite. In general 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 a multigrid method—details of this multigrid solver are presented later. A connection to the general continuous-time framework presented in Sec. 2.4 is provided in A.

3.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 $\Delta x=\Delta y=\Delta z=h$. The stress tensor $\boldsymbol{\sigma}$ and effective temperature $\chi$ are stored at cell centers and indexed by half-integers, while the velocity $\mathbf{v}$ is stored at cell corners and indexed by integers, as shown in Fig. 2(a).

Let $(\partial f/\partial x)_{i,j,k}$ denote the partial derivative of a field $f$ with respect to $x$ evaluated at grid point $(i,j,k)$. The staggered centered difference is

[TABLE]

Equation 31 averages four edge-centered centered differences surrounding the cell center and has a discretization error of size $\mathcal{O}(h^{2})$. The derivative at a cell corner is obtained by the replacement $(i,j,k)\rightarrow(i-\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2})$. The diffusive term appearing in the velocity update in Eq. 24 is computed via the standard centered difference formula,

[TABLE]

The advective derivatives in Eqs. 24 & 25 must be upwinded for stability; we use the second-order essentially non-oscillatory (ENO) scheme [91]. With $\left[f_{xx}\right]_{i,j,k}$ denoting the second derivative with respect to $x$ of the field $f$ at grid point $(i,j,k)$ computed using Eq. 32, the ENO derivative is defined in the $x$ direction as

[TABLE]

Equation 33 uses the curvature of $f$ to switch between an upwinded three-point derivative and a centered difference. Versions of Eqs. 31, 32, & 33 in the $y$ and $z$ coordinates are obtained analogously.

To solve Eq. 30 via numerical linear algebra, the spatial derivatives must first be discretized using finite differences. In addition to the finite differences discussed above, Eq. 30 also contains mixed partial derivatives. The $xy$-derivative is computed numerically as

[TABLE]

with analogous expressions for other mixed partials.

3.5 Parallelization via MPI and domain decomposition

We solve Eqs. 24 & 25 in parallel with a custom C++ code that uses the MPI library for parallelization [92]. We use a global grid comprised of $Q\times M\times N$ grid cells. Each grid cell is labeled by its lower corner index $(i,j,k)$ as shown in Fig. 2(a), with indices in the ranges $i\in\{0,1,\ldots,Q-1\}$, $j\in\{0,1,\ldots,M-1\}$, $k\in\{0,1,\ldots,N-1\}$. The global grid is split across processors using a standard Cartesian decomposition into $P_{x}\times P_{y}\times P_{z}$ subdomains. Specifically, each processor is indexed as $(I,J,K)$ with $I\in\{0,\ldots,P_{x}-1\}$, $J\in\{0,\ldots,P_{y}-1\}$, and $K\in\{0,\ldots,P_{z}-1\}$. Processor $(I,J,K)$ is responsible for the grid cells

[TABLE]

Hence the subdomain assigned to each processor is identical up to one grid point in each direction.

The finite difference stencils in Eqs. 31 & 32 require data from adjacent gridpoints, and the ENO derivative in Eq. 33 requires data from at most two grid points away. On grid points within two points of a subdomain boundary, the derivative calculation can therefore require inaccessible data in a distributed memory setting. To handle this, we pad each processor subdomain with ghost regions of width two. A ghost region is a cubical shell of non-physical grid points whose field values are filled with data from adjacent processors, so that each subdomain can locally access values computed and stored in adjacent subdomains (Fig. 2(b)). At the simulation boundaries, ghost regions are used to enforce boundary conditions. For periodic boundary conditions, the ghost regions wrap around to processor subdomains on the other side of the simulation. For non-periodic conditions, the ghost grid points store values necessary to enforce the given boundary condition: for a free boundary, the fields are linearly extrapolated into the ghost region. For a Dirichlet condition, the physical field values are reflected over the boundary into the ghost region.

At the start of each timestep, each processor communicates with the 26 nearby processors sharing six faces, twelve edges, and eight corner regions. Each processor sends data to nearby processors and receives the data it requires from the same nearby processors via non-blocking communication. Each processor waits to ensure it has received all face data before updating faces, edge data before updating edges, and corner data before updating corners. This ensures that components of each ghost region can be updated when the data is available without waiting for all data to arrive. The total cost of parallel communication scales with the surface area of a processor subdomain. Usually, the total number of processors $P_{\textrm{tot}}$ is chosen to be power of two, and the precise subdomain decomposition is determined by considering all triplets $(P_{x},P_{y},P_{z})$ that satisfy $P_{x}P_{y}P_{z}=P_{\textrm{tot}}$, and selecting the one that minimizes the surface area between the subdomains.

3.6 Performing the projection step

We solve Eq. 30 for the velocity using a custom parallel implementation of the geometric multigrid method, a linear system solver that is particularly suited to elliptic problems that take place on a physical grid [27]. Let $G_{0}$ be the original grid, and let $\mathbf{A}_{0}\mathbf{x}_{0}=\mathbf{b}_{0}$ be the linear system to solve on this grid. In the multigrid method, a hierarchy of progressively coarser grids $G_{1},G_{2},\ldots,G_{g}$ is introduced. In our implementation, if $G_{k}$ has resolution $Q_{k}\times M_{k}\times N_{k}$, then $G_{k+1}$ has resolution $\lceil Q_{k}/2\rceil,\lceil M_{k}/2\rceil,\lceil N_{k}/2\rceil$. Interpolation operators $\mathbf{T}_{k}:G_{k}\to G_{k-1}$ are introduced based on linear interpolation, and restriction operators $\mathbf{R}_{k}:G_{k}\to G_{k+1}$ are introduced based on local averaging. Both $\mathbf{T}_{k}$ and $\mathbf{R}_{k}$ can be represented as rectangular matrices, and in our implementation $\mathbf{R}_{k-1}=\mathbf{T}_{k}^{\mathsf{T}}$—this condition is not necessary for a practical implementation, but is useful in some convergence proofs [28].

Our multigrid implementation uses the standard V-cycle [27, 28] with two pre-smoothing steps and two post-smoothing steps. On $G_{0}$ the grid is decomposed among the processors in the same way as the simulation fields. The smoothing steps are performed using the Gauss–Seidel method on each processor, with the ghost regions being synchronized after each step. This requires building a representation of the linear system on each grid, which we do via recursive matrix multiplication [93, 94],

[TABLE]

The implementation works with periodic and non-periodic boundary conditions, and arbitrary grid dimensions. As the grids are coarsened, the amount of work on each grid is rapidly reduced, to the point where it is no longer effective for all processors to share the work. The implementation therefore has the ability to amalgamate the coarser problem onto a smaller set of processors, with the rest remaining idle.

The multigrid implementation uses C++ templates, so that the linear system can be compiled to work with an arbitrary data type. For the current problem, $\mathbf{b}_{0}$ is given by the source term $-\Delta t\nabla\cdot\boldsymbol{\sigma}$ and $\mathbf{x}_{0}$ contains values of $\mathbf{v}^{n+1}$ across the entire grid. Hence, we compile the multigrid library where the elements of $\mathbf{b}_{0}$ and $\mathbf{x}_{0}$ are 3-vectors, and the block elements of $\mathbf{A}_{0}$ are $3\times 3$ symmetric matrices. The matrix $\mathbf{A}_{0}$ is sparse, and a grid point $(i,j,k)$ is only coupled to the 27 grid points in the $3\times 3\times 3$ surrounding cube of grid points given by coordinates $(i+\{-1,0,1\},j+\{-1,0,1\},k+\{-1,0,1\})$ in our discretization scheme.

4 Shearing between two parallel plates

In the following sections, we consider several material samples being sheared between two parallel plates. This example is experimentally relevant, has simple boundary conditions, demonstrates complex shear banding dynamics [76, 74, 73, 30, 77, 49, 48, 75], and has been studied previously in two dimensions [13]. It represents a natural physical scenario to compare three-dimensional results to two-dimensional results, compare simulation data to experiments, and to quantitatively compare the explicit and quasi-static methods.

In all simulations, we consider a domain periodic in the $x$ and $y$ directions with shear velocity applied on the top and bottom boundaries in $z$. The domain occupies $-L\leq x<L,-L\leq y<L$, and $-\gamma L\leq z\leq\gamma L$ with $\gamma=\frac{1}{2}$ and $L=1\text{~{}cm}$. A natural unit of time is given as $t_{s}=L/c_{s}$ where $c_{s}=\sqrt{\mu/\rho}$ is the material shear wave speed, and we measure time in this scale. The boundary conditions are given by

[TABLE]

where the function $U(t)$ is given by

[TABLE]

The ramp-up in the function $U(t)$ prevents a large deformation rate near the boundary that would be present with $U(t)=U_{B}$ immediately at $t=0$. The elasticity and plasticity parameters are defined in Table 1. From these values, the natural timescale is $t_{s}=4.05~{}\text{\textmu{}s}$.

A diagram of the global three-dimensional grid and the ghost regions at simulation boundaries used for implementing the boundary conditions is shown in Fig. 3. The cell-cornered grid points run according to $i\in\{0,\ldots,Q-1\}$, $j\in\{0,\ldots,M-1\}$, and $k\in\{0,\ldots,N\}$; because the grid is non-periodic in $z$ there is an extra grid point in this direction. The velocities at grid 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\in\{\frac{1}{2},\frac{3}{2},\ldots Q-\frac{1}{2}\}$, $j\in\{\frac{1}{2},\frac{3}{2},\ldots M-\frac{1}{2}\}$, and $k\in\{\frac{1}{2},\frac{3}{2},\ldots N-\frac{1}{2}\}$. Ghost layers of cell-centered grid points are at $(i,j,-\frac{1}{2})$, $(i,j,-\frac{3}{2})$, $(i,j,N+\frac{1}{2})$, and $(i,j,N+\frac{3}{2})$ with equivalent expressions in the other dimensions and for cell-cornered grid points. The values of $\boldsymbol{\sigma}$ and $\chi$ in the ghost layers are linearly extrapolated from the two nearest layers to ensure that these fields remain free on the boundary. At the simulation boundaries in the $x$ and $y$ directions, ghost points outside the simulation domain are filled with values that wrap around.

Following the introduction in Sec. 3.1, the effective temperature $\chi$ is a continuum-scale variable that encodes the density of STZs, and hence the details of the material’s microscopic structure. Its initial condition determines the evolution of plastic deformation within the material.

4.1 Qualitative comparison between explicit and quasi-static methods

We now demonstrate the qualitative equivalence between results computed with the explicit and quasi-static methods. We consider an initial condition corresponding to a finite cylindrical inclusion

[TABLE]

Initially the cylindrical inclusion is slightly more amenable to plastic deformation, and hence we expect to see a shear band nucleate from it. To visualize the effective temperature field in three dimensions, we use a custom opacity function defined as

[TABLE]

Equation 40 sets the opacity of a grid point based on the value of $\chi(\mathbf{x})$. The parameters $a$ and $\eta$ are chosen on a case-by-case basis to reveal the most interesting features111Ideally, we would like to use the same opacity parameters for all plots. However, due to significant variations in the ranges of the $\chi$ fields and in their spatial structures, we found it was necessary to set the parameters on an individual basis. We note, however, that the color scale is absolute across all simulations.. The initial condition is depicted in Fig. 4. The grid is of size $64\times 64\times 32$ to accommodate the limitations of the explicit simulation method, corresponding to a grid spacing $h=\frac{L}{32}$. A viscous stress constant of $\kappa^{\prime}=4.8\frac{L}{t_{s}}$ is used in Eq. 24. The timestep for the explicit method is $\Delta t_{\text{e}}=\frac{h^{2}t_{s}}{20L^{2}}$.

A typical applied shear velocity that is comparable to a realistic loading rate in a laboratory experiment is $U_{b}=10^{-7}L/t_{s}$ [76]. With this velocity, running an explicit simulation is prohibitively expensive due to the CFL condition. To ensure that significant plastic deformation occurs on timescales reachable by the explicit method, a scaling parameter $\zeta$ is introduced. The molecular vibration timescale $\tau_{0}$ is rescaled to $\tau_{0}\zeta^{-1}$ and the applied shear velocity is inversely rescaled to $U_{B}=10^{-7}L/t_{s}\zeta$. The simulation is conducted until a final time of $t_{f}=2\times 10^{6}t_{s}/\zeta$. As $\zeta$ approaches zero, the quasi-static limit of Eqs. 1 & 2 is formally approached. We therefore expect greater agreement for smaller values of $\zeta$. Due to the appearance of $\frac{1}{\tau_{0}}$ in Eq. 20, the introduction of $\zeta$ has the effect of linearly scaling the magnitude of plastic deformation by a factor of $\zeta$. A quasi-static timestep of $\Delta t_{\text{qs}}=200t_{s}/\zeta$ is used.

In Fig. 5, we show three snapshots of the effective temperature field with $\zeta=10^{4}$ from each of the two simulation methods, at $t=50t_{s}$, $t=75t_{s}$, and $t=100t_{s}$ respectively. The explicit 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=50t_{s}$, a shear band begins to emerge, nucleating outwards from the center of the simulation. A thin region of higher $\chi$ is visible in the center of the band. By $t=75t_{s}$, the shear band has fully formed and spans the system. At $t=100t_{s}$, the band grows stronger and $\chi$ continues to increase.

Figure 6 shows cross-sections in $z$ for fixed $x=0$ and $y=0$ of $\|\boldsymbol{\sigma}_{0}\|_{\text{qs}}-\|\boldsymbol{\sigma}_{0}\|_{\text{e}}$ for several time points before the onset of plastic deformation, highlighting some differences between the two methods. The explicit simulation exhibits oscillations due to elastic waves propagating through the medium. Because the quasi-static method does not resolve these elastic waves, the oscillations are apparent in the deviatoric stress differences. When plastic deformation sets in, plasticity-induced damping removes the elastic waves and the agreement improves.

4.2 Quantitative comparison between explicit and quasi-static methods

Having demonstrated the qualitative agreement between the two simulation methodologies for $\zeta=10^{4}$ in the previous section, we now examine convergence as $\zeta$ is decreased. The same simulation geometry, boundary conditions, and initial conditions in the effective temperature field are used here as in the previous section. To quantitatively compute the agreement between the explicit and quasi-static methods, we define a norm on simulation fields $\mathbf{f}$,

[TABLE]

The integral in Eq. 41 runs over the entire simulation domain and is computed numerically via the trapezoid rule. The appearance of $|\cdot|$ in Eq. 41 is taken to be the Euclidean norm for vectors, absolute value for scalars, and the Frobenius norm for matrices.

Equation 41 is evaluated for $\chi_{\text{qs}}-\chi_{\text{e}}$, $\boldsymbol{\sigma}_{\text{qs}}-\boldsymbol{\sigma}_{\text{e}}$, and $\mathbf{v}_{\text{qs}}-\mathbf{v}_{\text{e}}$ at intervals of $0.02t_{s}$ in pairs of explicit and quasi-static simulations with $\zeta=10^{4},5\times 10^{3},2.5\times 10^{3}$ and $1.25\times 10^{3}$. In each case, the norm is non-dimensionalized using the quantities $\chi_{\infty}$, $s_{Y}$, and $U_{B}$ to ensure the norm values for quantities of different units are of comparable magnitude. For each simulation, the explicit timestep is $\Delta t_{\text{e}}=\frac{t_{s}h^{2}}{20L^{2}}$ and the quasi-static timestep is $\Delta t_{\text{qs}}=\frac{100t_{s}}{\zeta}$.

Plots of all three norm values are shown as a function of time in Fig. 7(a) for the value of $\zeta=1.25\times 10^{3}$. Oscillations due to elastic waves are visible in all simulation fields until around $t\approx 12t_{s}$ when the yield stress is reached. After the onset of plastic deformation, the norm in effective temperature increases steadily, most rapidly during the period of shear band nucleation from $t\approx 12t_{s}$ to $t\approx 25t_{s}$. The disagreements in $\boldsymbol{\sigma}$ and $\mathbf{v}$ decrease during the elastic region, and steadily increase after plastic deformation begins.

The remaining three panels in Fig. 7 show the quantitative comparisons as a function of time for values of $\zeta=10^{4},5\times 10^{3},2.5\times 10^{3}$, and $1.25\times 10^{3}$ for $\mathbf{v}$, $\chi$, and $\boldsymbol{\sigma}$. In all plots, better agreement with smaller $\zeta$ is observed during the elastic regime and during the onset of plasticity while $t\leq 12t_{s}$. After shear band nucleation from $12t_{s}\leq t\leq 25t_{s}$, all values of $\zeta$ have roughly equal error magnitudes, with slightly greater agreement for higher values of $\zeta$. This is consistent with previous comparisons in two dimensions, where the dominant factor governing the disagreement between the two simulation methods was shown to be due to differences in the discretization rather than the value of $\zeta$ itself [13].

A method to reduce the differences in discretization is to increase the background $\chi$ field. With higher values of background $\chi$, finer-scale features in the shear banding dynamics are less 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 timestep between the two methods will be less significant. Snapshots of the effective temperature field are shown in Fig. 8 at $t=10^{6}t_{s}$ for background $\chi$ field values of $\chi_{\text{bg}}=600\text{~{}K},650\text{~{}K}$, and $700\text{~{}K}$. Figure 9 confirms that the differences between the two types of simulation decreases as $\chi_{\text{bg}}$ increases.

4.3 Parallel scaling analysis

In this subsection, we present a parallel scaling analysis of our quasi-static algorithm implementation. All simulations are run on an Ubuntu Linux computer with dual 16-core 2.10 GHz Intel Xeon E5-2683 v4 processors. We consider two measures of parallel scaling. First, we vary the number of processors $p$ while keeping $N_{x}N_{y}N_{z}/p$ approximately fixed as a test of weak parallel scaling. We denote $N=N_{x}=N_{y}=2N_{z}$ so that the total number of grid points is $N^{3}/2$. For this weak scaling experiment, we consider an initial condition in the effective temperature field corresponding to a non-commensurate helix, given by

[TABLE]

We consider the values $N=80,102,128,164,204,256$ and matching values $p=1,2,4,8,16,32$. Each experiment uses a boundary shear velocity of $U_{b}=10^{-7}\frac{L}{t_{s}}$, a diffusion lengthscale of $l=\frac{3}{2}h$, a quasi-static timestep of $\Delta t=200t_{s}$, and is simulated to a final time $t_{f}=10^{6}t_{s}$. The shear banding dynamics for $N=256$ are shown in Fig. 10. In Fig. 10(a), the initial condition is shown. In Fig. 10(b), the background effective temperature field begins to increase, but shear banding has not yet begun. In Fig. 10(c), clear system-spanning shear bands have begun to nucleate off the top and bottom of the helix. In Figs. 10(d)–(f), the shear bands continue to become more prominent and thicken.

Experimental results are shown in Table 2. The weak scaling ratio is computed using the formula $w=\frac{rt_{1}}{t_{p}}$ where $t_{1}$ is the time for the simulation with one processor, $t_{p}$ is the time for the simulation with $p$ processors, and $r=N^{3}_{p}/(pN_{1}^{3})$ is a small adjustment factor for deviations from perfectly equal work. The total time increases with grid size at a significantly reduced rate when compared to serial computation, and the weak scaling ratios for $p\neq 32$ indicate good weak scaling performance. The value of the weak scaling ratios at $p=32$ may indicate an issue such as limitations of communication buffers. These results hold for both the total time and the V-cycle computation, with the V-cycle computations showing mildly worse weak scaling than the overall computation. The number of required V-cycles increases slightly as the grid size is increased.

To quantitatively determine the parallel scaling efficiency, we assume that the total work required can be written as $W=k\times p$ where $k$ is a constant. We further assume that the total time taken can be computed as $T=W/(pQE(p))$ where $Q$ is the total power per processor and $E(p)$ describes the parallel efficiency. With $E(p)\propto p^{-\alpha}$ for some exponent $\alpha$, the logarithm of the total time taken is linear in $\log(p)$ with coefficient $\alpha$, $\log T=A+\alpha\log p$ where $A$ is a constant, so that $\alpha$ may be found via linear regression. We have computed the exponent $\alpha$ for the total time taken and the time per V-cycle, as shown in Fig. 11. The $\alpha$ value in each case was found to be $\alpha=0.2258$ and $\alpha=0.2092$.

We also considered testing strong scaling performance of our implementation. However, fixing the size of the grid and varying the number of processors is infeasible in this setting. If the simulation is small enough that it can be handled in a reasonable amount of time, going to higher processor counts will not be beneficial due to communication limitations. If the simulation is large enough that higher processor counts such as $p=32$ are useful, the simulation will take too long to run with few processors.

As an alternative measure of parallel efficiency, we consider fixing $p=32$ and varying the number of grid points $N^{3}/2$. In this experiment, we use the same simulation parameters as for the test of weak scaling, but now with a random (drawn from a normal distribution) initial condition for the effective temperature with mean $550\text{~{}K}$ and standard deviation $15\text{~{}K}$; more details on the initialization procedure are provided in Sec. 4.4.2. The results listed in Table 3 demonstrate that the number of V-cycles per simulation again increases with $N$. The total time per simulation and time per V-cycle increases greater than linearly. Fitting $\log T=B+\beta\log N$ in Fig. 11, we find that $\beta=2.395$ for the total simulation time and $\beta=2.275$ for the average time per V-cycle. Over the range of sizes considered, our implementation has more favorable scaling than the predicted $\mathcal{O}(N^{3})$ of multigrid. This is because the communication becomes comparatively more efficient for larger grid sizes, since the ratio of ghost points to total grid points tends to zero as $N$ increases.

4.4 Gaussian defects

We now turn to simulating realistic physical timescales with the quasi-static method, where we now set the scaling parameter $\zeta=1$. We first consider the nucleation of shear bands from localized imperfections of higher $\chi$. Physically, these imperfections describe defects within the material structure that may be particularly susceptible to plastic deformation [95]. To begin, we consider a single defect, corresponding to an initial $\chi$ field of the form

[TABLE]

The simulation is performed on a grid of size $256\times 256\times 128$, corresponding to a grid spacing of $h=L/128$. The length scale $l$ appearing in Eq. 23 is fixed at $h$ and sets the width of the shear bands. The boundary velocity is set to a value of $U_{b}=10^{-7}L/t_{s}$, and the simulation is conducted to a final time of $t_{f}=10^{6}t_{s}$, using a quasi-static timestep of $\Delta t=200t_{s}$. For three-dimensional visualization, we use the opacity function from Eq. 40. The simulation takes $7.52$ total hours using 32 processes on an Ubuntu Linux computer with dual 10-core Intel Xeon E5-2630 v4 processors. The total time spent in V-cycle computations is $5.852$ hours. The total number of V-cycles is $17586$.

Snapshots of the effective temperature field at various time points are shown in Fig. 12. The initial condition is shown in Fig. 12(a). At $t=2\times 10^{5}t_{s}$ in Fig. 12(b), the defect has started to expand. By $t=4\times 10^{5}t_{s}$ in Fig. 12(c), a shear band begins to nucleate, indicated by a quadrupolar structure emanating from the defect. The background $\chi$ field also begins to increase, as demonstrated by the presence of the transparent light blue background. By $t=6\times 10^{5}t_{s}$ in Fig. 12(d), a distinct system-spanning band has become clear. The band displays no curvature in either of the $x$ or $y$ directions. By $t=8\times 10^{5}t_{s}$ in Fig. 12(e), a prominent band has formed. $t=10^{6}t_{s}$ in Fig. 12(f) is similar, with a thicker and stronger band.

We now introduce a second defect to highlight some three-dimensional characteristics of shear banding. We expect that the relative size and displacement between the defects will determine the dynamics, with the possibility of forming a single shear band that connects the two. The initial effective temperature field is

[TABLE]

Two cases of Eq. 44 are considered. First, we take $\mathbf{X}_{1}=(-0.5,-0.5,0.35)$ and $\mathbf{X}_{2}=(-0.5,0.5,0.25)$, corresponding to two defects symmetric about the $y=0$ plane with the same $x$ coordinate, a slight offset in $z$, and different sizes. The results for this case using the same simulation parameters as Fig. 12 are shown in Fig. 13. Second, we take $\mathbf{X}_{1}=(-0.5,-0.5,0.35)$ and $\mathbf{X}_{2}=(0.5,-0.5,0.25)$; this is the same as the previous case, but with the roles of $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. 14. The initial configurations are shown in Figs. 13(a) and 14(a). The first arrangement takes $5.61$ total hours using 32 processes on a Ubuntu Linux computer with dual 10-core 2.20 GHz Intel Xeon E5-2630 v4 processors. The total time spent in multigrid V-cycles is $4.25$ hours. The total number of V-cycles is $17782$. The second arrangement takes $5.27$ hours using 32 processes on an Ubuntu Linux computer with dual 10-core 2.20 GHz Intel Xeon Silver 4114 processors. The total time spent in multigrid V-cycles is $3.98$ hours. The total number of V-cycles is $17647$.

There is an interesting contrast between the time sequences displayed in Figs. 13 and 14. Much like in the single defect simulations, $t=2\times 10^{5}t_{s}$ in Figs. 13(b) and 14(b) displays expansion of the defects, and $t=4\times 10^{5}t_{s}$ in Figs. 13(c) and 14(c) shows the initiation of shear band nucleation. In Fig. 13(d), we see the formation of a single curved band connecting the two defects, while in Fig. 14(d), the band is flat. Figs. 13(e), 13(f), 14(e), and 14(f) make this more clear as the band becomes more defined. The curvature seen in Fig. 13 is in the direction orthogonal to shear.

The dependence of band curvature on the relative orientation of the two defects can be best understood in terms of the qualitative structure of Fig. 12(c). There is a substantial extension of elevated $\chi$ along the $x$ direction (parallel to shear, in-plane), a small extension along the $y$ direction (orthogonal to shear, in-plane), and a moderate extension along the $z$ direction (orthogonal to shear, out-of-plane). In Fig. 13, the defects are offset in $y$ and $z$. Because the $\chi$ field is stronger in $z$ than in $y$, this relative placement of the defects can accomodate curvature along the $y$ direction. On the other hand, in Fig. 14, the defects are offset in $x$ and $z$. The strength of the $\chi$ field extension in the $x$ direction is great enough that the flat, horizontal regions of the two forming bands reach each other. The two bands thus join into one fatter flat band.

Taken together, Figs. 12, 13, & 14 provide insight into the structure and nucleation of shear bands from localized material defects. They help understand experimentally observed band curvature and raise the possibility that the placement and orientation of microscopic material properties can influence the qualitative structure of macroscopic shear bands. Finally, they provide intuition for more complex initial conditions, such as the random initializations considered later in this work, as a superposition of many defects.

4.4.1 Circular defects

We now turn to a set of more complex initial conditions in the effective temperature field. Results for initial conditions corresponding to a circular region of elevated $\chi$ parallel and perpendicular to the direction of shear are shown in Figs. 15 and 16 respectively. The initial conditions are given by

[TABLE]

and

[TABLE]

representing circles in the $yz$ and $xz$ planes respectively. Simulations are carried out using the same simulation geometry, discretization, quasi-static timestep, and boundary velocity as in the previous section. The initial conditions in Eqs. 45 and 46 are displayed in Figs. 15(a), 16(a) respectively. The first arrangement takes $5.66$ total hours using 32 processes on a Ubuntu Linux computer with dual 14-core 1.70 GHz Intel Xeon E5-2650L v4 processors. $4.30$ hours are spent in multigrid V-cycles. There are $16980$ total V-cycles. The second arrangement takes $5.46$ total hours using 32 processes on an Ubuntu Linux computer with dual 8-core 2.40 GHz Intel Xeon E5-2630 v3 processors. $4.22$ hours are spent in multigrid V-cycles. There are $17153$ total V-cycles.

By $t=2\times 10^{5}t_{s}$ in Figs. 15(b) and 16(b), little has changed, though the circles have expanded slightly. At $t=4\times 10^{5}t_{s}$ in Figs. 15(c) and 16(c), differences due to the orientation of the circles become clear. The circle oriented perpendicular to shear closes vertically into a disk. The circle oriented along shear exhibits signatures of shear band nucleation at four equally-spaced points. At $t=6\times 10^{5}t_{s}$ in Fig. 15(d), the disk has expanded and has developed two diffuse shear bands. In Fig. 16(d), there are two thinner, sharper, more well-separated bands forming off the top and bottom of the circle. By $t=8\times 10^{5}t_{s}$ in Figs. 15(e) and 16(e), these differences have become even more prominent. The bands are seen to have a curved structure in the $y$ direction in Fig. 15(e) which is not as clear in Fig. 16(e). The interim region is of higher $\chi$ than for the bands seen in Fig. 16(e), where the circular defect is still mostly visible. These features continue to develop into the final pane at $t=10^{6}t_{s}$. Taken together, Figs. 15 and 16 demonstrate another example of the dependence of shear banding structure and dynamics on the orientation of initial conditions in the $\chi$ field with respect to shear.

4.4.2 A randomly fluctuating effective temperature field

In this section, we consider the case of a randomly distributed initial effective temperature field. The initial conditions presented in the previous sections provide insight into the dynamics of shear banding, but it is unlikely that they have exact physical correspondences. The STZ theory postulates that STZs are randomly distributed throughout the material, and a random initial condition in $\chi$ is most faithful to this fundamental assumption [96]. Random initial conditions are thus expected to shed the most light on the structure of shear bands observed in experiments. The randomly fluctuating $\chi$ field leads to the formation of multiple shear bands, and may enable the study of shear band interactions in the STZ model [97].

We first populate the grid and a shell of ghost points with random variables $\chi_{\zeta}(\mathbf{x})$ using the Box–Muller algorithm. With $\mu_{\chi}$ and $\sigma_{\chi}$ respectively denoting the desired mean and standard deviation, we perform the convolution

[TABLE]

where $V$ denotes the set of grid points and $V^{\prime}$ denotes the set of grid points with the addition of the ghost points. Equation 47 ensures that the effective temperature value at each point is normally distributed with mean $\mu_{\chi}$ and standard deviation $\sigma_{\chi}$. In practice, the sums in Eq. 47 are performed with a cutoff length scale specified as a multiplicative factor of the convolution length scale $l_{c}$, and the number of ghost points in $V^{\prime}$ is set by the choice of cutoff length scale. Results for a random initialization with $\mu_{\chi}=550\text{~{}K}$, $\sigma_{\chi}=15\text{~{}K}$, $l_{c}=10h$ and a cutoff factor of $5$ (leading to $50$ ghost points in each direction for the convolution) are shown in Fig. 17. The grid is of size $768\times 768\times 384$. The simulation geometry, quasi-static timestep, and boundary velocity are the same as in previous sections. The initial conditions are shown in Fig. 17(a). The simulation takes 105.970 total hours using 32 processes on an Ubuntu Linux computer with dual 16-core Intel Xeon E5-2683 v4 processors. The total time spent in V-cycles is 75.90 hours. The total number of V-cycles is 19689.

By $t=2\times 10^{5}t_{s}$ in Fig. 17(b), the effective temperature has increased somewhat uniformly across the grid. At $t=4\times 10^{5}t_{s}$ in Fig. 17(c), both horizontal and vertical shear bands begin to nucleate throughout the simulation. Slightly later at $t=6\times 10^{5}t_{s}$ in Fig. 17(d), a multitude of thin, system-spanning horizontal bands connected by vertical bands have begun to emerge. Curvature is present in the horizontal bands both parallel and orthogonal to the direction of shear. The shear bands become increasingly prominent and grow in number by $t=8\times 10^{5}t_{s}$ and $t=10^{6}t_{s}$ in Figs. 17(e) and (f) respectively, where crossing and branching patterns in the many bands are observed.

5 Friction welding

In this section, we investigate the shear banding dynamics for a set of boundary conditions inspired by friction welding, a process of physical interest and engineering relevance [98, 99, 100]. A typical friction welding example consists of pushing a rapidly spinning cylinder into a block of material. The friction generated at the block–cylinder interface causes the material to heat up and melt, forming a weld. Inspired by this scenario, we examine how applying spinning-disc boundary conditions to a block of material can generate deformation. Let $r^{2}=\left(\frac{x}{L}\right)^{2}+\left(\frac{y}{L}\right)^{2}$. We impose a boundary condition corresponding to

[TABLE]

[TABLE]

on the top boundary. $\lambda_{w}$ is a parameter that sets the smoothness of the boundary condition to avoid stress concentration at the boundary of the rotating region; in the limit $\lambda_{w}\rightarrow\infty$, Eqs. 48&49 describe a rotation only within the disc $r^{2}<R^{2}$. In all simulations, we set $\lambda_{w}=16\log(10)/(1-R^{2})^{2}$ and $R=\frac{1}{3}$. This choice of $\lambda_{w}$ ensures that one quarter of the way between the boundary of the spinning disc $r^{2}<R^{2}$ and the edge of the simulation, the angular velocity is $10\%$ of its value inside the disc.

To promote shear banding with this boundary condition, we set $\gamma=\frac{1}{4}$, corresponding to a more slender domain than considered in previous sections. All simulations are conducted on a $768\times 768\times 384$ cell grid with a random initial effective temperature distribution as described in the preceding section. We choose values $\mu_{\chi}=450$ and $\sigma_{\chi}=15$. The convolution lengthscale and cutoff factor are set to $5h$ and $5$ respectively. We simulate until a final time of $t_{f}=10^{8}t_{s}$ and use a quasi-static timestep of $\Delta t=2\times 10^{4}t_{s}$. We choose a value of $\omega=\frac{2\pi}{t_{f}}$, corresponding to a single rotation over the course of the simulation. We found that lowering the angular velocity by increasing the simulation duration promoted interesting shear banding dynamics, and we increased the timestep accordingly.

Three-dimensional visualizations of the effective temperature distribution are shown in Fig. 18. In these figures, it is clear that the plastic deformation is predominantly concentrated on the top boundary, and but propagates a short distance below the disc into the bulk. To better understand the structure of this leakage, we consider $xz$ cross sections of the $\chi$ field for fixed $y$ in Fig. 19 at $t=5\times 10^{7}t_{s}$. These cross sections are shown at eight roughly equally spaced points on one side of the $y=0$ plane; the dynamics are symmetric up to discrepancies in the random initialization. Here, we see a parabolic region of elevated $\chi$ that begins to pucker as the center of the disc is approached.

Given the structure of the $\chi$ field in the cross sections in Fig. 19, we may expect more interesting dynamics if both the top and bottom surfaces are welded. In Fig. 20, we show three-dimensional visualizations of the shear banding dynamics when Eqs. 48 and 49 are applied to both the top and bottom boundaries. Here, the deformation dynamics are much richer, and complex structure unlike that produced by simple shear can be seen. In Fig. 20(b), the discs of elevated $\chi$ on the top and bottom boundaries have both begun to propagate into the bulk. Even at this stage, the leakage is more complicated than in Fig. 18, and includes thin extensions of elevated $\chi$. In Fig. 20(c)–(f), the structures on the top and bottom begin to connect with each other, forming a hyperboloidal region of elevated $\chi$ spanning the system in $z$.

To probe the structure of this hyperboloid, we again consider two-dimensional cross sections in Fig. 21 at $t=5\times 10^{7}t_{s}$. The hyperboloid does not extend entirely through the bulk, but instead has a gap between the two regions of elevated $\chi$ created by the spinning discs. In the gap are thin, tendril-like shear bands that work to connect the two regions and close the gap.

Finally, we may ask how the structure will change if the top and bottom are anti-rotated with respect to each other. Three-dimensional visualizations for this setup are shown in Fig. 22. A similar hyperboloid region to Fig. 20 is seen, though it develops later in the simulation, and does not occupy as much of the disc $r^{2}<R^{2}$. This structure is made more clear in the contour plots in Fig. 23, where the connecting tendrils only form near the center of the disc and are thinner.

6 Conclusion

Expanding on prior two-dimensional work [13], we have developed a three-dimensional numerical method for simulating quasi-static elastoplastic materials by analogy with the projection method of Chorin for incompressible fluid dynamics [31, 32]. The method is particularly suitable for stiff materials with small elastic deformation and high elastic wave speeds, where plastic deformation is often quasi-static, and the hypo-elastoplastic assumption is 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 the quasi-static method has no such restriction.

The method was used to examine the properties of shear banding in a metallic glass with the STZ plasticity model, focusing on uniquely three-dimensional features such as curvature development, nucleation off of complex geometries, dependence of shear band structure on the orientation of features in the initial $\chi$ distribution, and high-resolution banding with a random initialization in $\chi$. Essential to these studies was the development of a high performance, parallelized, three-dimensional geometric multigrid solver, along with an MPI and C++-based implementation of the projection algorithm. We provided a weak parallel scaling analysis for our implementation and showed that it exhibits favorable weak scaling over a range of processor counts. The method, though applied to the study of BMGs in this work, is independent of the plasticity model.

The method is based on a correspondence between the variables $(\mathbf{v},p)$ in incompressible fluid dynamics and the variables $(\boldsymbol{\sigma},\mathbf{v})$ in quasi-static hypo-elastoplasticity. As part of the development of the method, we introduced an auxiliary vector field $\boldsymbol{\Phi}$, which plays a role analogous to the auxiliary scalar field $\phi$ used in projection algorithms for fluid dynamics. The choice of $\phi$, along with careful consideration of its boundary conditions, leads to higher-order projection and gauge methods in fluid dynamics. The analogy between the two auxiliary fields is reminiscent of the analogy between the physical fields that led to the development of the quasi-static method, and it suggests that similar approaches may generalize to quasi-static hypo-elastoplasticity [71, 68, 69, 70]. Gauge methods couple naturally with discontinuous Galerkin discretization (dG) methods, and the theoretical developments presented here thus pave the way for the possibility of similar dG methods in quasi-static hypo-elastoplasticity. Such methods may help resolve the fine-scale features of instabilities like shear bands, and the development of these higher-order methods is an interesting avenue for 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.4 as follows. By comparison of Eqs. 17 and 30, we can identify $\boldsymbol{\Phi}=\Delta t\ \mathbf{v}$. Equation 16 then says that

[TABLE]

Recall that $\mathbf{q}$ is chosen to be the best available approximation to $\mathbf{v}^{n+1}$, and note that by symmetry of $\mathbf{C}$,

[TABLE]

Equation 50 thus says the following: $\mathbf{C}:\mathbf{D}^{n+1}$ is given by the best available guess before the solve for $\mathbf{v}^{n+1}$ - $\mathbf{C}:\nabla\mathbf{q}$ - plus an $\mathcal{O}(\Delta t)$ correction constructed via a first-order Taylor expansion 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 Science Foundation under Grant Nos. DMR-1409560 and DMS-1753203. N. M. Boffi was supported by a Department of Energy Computational Science Graduate Fellowship under Grant No. DE-FG02-97ER25308. C. H. Rycroft was partially supported by the Applied Mathematics Program of the U.S. DOE Office of Advanced Scientific Computing Research under contract number DE-AC02-05CH11231.