Coordinate transformation methodology for simulating quasi-static elastoplastic solids

TL;DR

This paper introduces a coordinate transformation method for simulating quasi-static elastoplastic solids, enabling precise comparison between molecular dynamics and continuum mechanics through a fixed grid framework.

Contribution

It develops a novel projection-based simulation framework with coordinate transformations to accurately model quasistatic elastoplastic deformation.

Findings

Effective simulation of shear band evolution in metallic glasses.

Comparison of shear and pure shear conditions on deformation.

Validation of the coordinate transformation approach.

Abstract

Molecular dynamics simulations frequently employ periodic boundary conditions where the positions of the periodic images are manipulated in order to apply deformation to the material sample. For example, Lees-Edwards conditions use moving periodic images to apply simple shear. Here, we examine the problem of precisely comparing this type of simulation to continuum solid mechanics. We employ a hypoelastoplastic mechanical model, and develop a projection method to enforce quasistatic equilibrium. We introduce a simulation framework that uses a fixed Cartesian computational grid on a reference domain, and which imposes deformation via a time-dependent coordinate transformation to the physical domain. As a test case for our method, we consider the evolution of shear bands in a bulk metallic glass using the shear transformation zone theory of amorphous plasticity. We examine the growth of shear bands in simple shear and pure shear conditions as a function of the initial preparation of the bulk metallic glass.

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 taken, the total amount of time spent in multigrid V-cycles, the total number of multigrid V-cycles, and the average time per V-cycle for the two simulation approaches. “T” specifies the transformed simulation and “NT” the non-transforemd simulation. The transformed method takes longer than the non-transformed method in general due to an increased number of multigrid V-cycles required to achieve convergence. The average time spent per V-cycle is roughly the same in the two approaches. Each simulation uses 32 32 32 processes.

	$N=64$		$N=96$		$N=128$		$N=160$		$N=192$		$N=256$
	T	NT	T	NT	T	NT	T	NT	T	NT	T	NT
Total time (hours)	0.0633	0.0343	0.5623	0.2863	2.4283	1.1981	3.3058	1.7890	8.5285	4.5787	33.3239	20.0178
V-cycle time (hours)	0.0452	0.0280	0.4130	0.2434	1.7845	1.0365	2.1976	1.3251	5.7242	3.4398	21.2268	15.2432
# of V-cycles	5544	3603	12481	8106	22181	14402	34658	22502	49913	32404	73164	57600
Time/V-cycle (seconds)	0.0294	0.0280	0.1191	0.1081	0.2896	0.2591	0.2282	0.2120	0.4129	0.3810	1.0444	0.9527

Table 3. Table 3: Data describing the total time, total time spent in multigrid V-cycles, total number of multigrid V-cycles, average time spent per multigrid V-cycle, and processor details for each randomly initialized simulation. This data applies to the randomly initialized simulations with non-periodic boundary conditions and simple shear deformation. The number of required multigrid V-cycles decreases as the background χ 𝜒 \chi field increases, likely due to more homogeneous dynamics. Each simulation uses 32 32 32 processes.

	${\mu }_{\chi }=450\text{ K}$	${\mu }_{\chi }=500\text{ K}$	${\mu }_{\chi }=525\text{ K}$	${\mu }_{\chi }=550\text{ K}$	${\mu }_{\chi }=575\text{ K}$	${\mu }_{\chi }=600\text{ K}$
Total time (hours)	95.2948	89.7239	76.9704	82.7853	71.7865	69.1593
V-cycle time (hours)	65.7853	60.5694	48.6470	53.4612	41.7683	40.4283
# of V-cycles	34846	30663	26628	24991	22649	20735
Time/V-cycle (seconds)	6.7964	7.1111	6.5749	7.7012	6.6390	7.0195
Processor details	Dual 10-core 2.20 GHz Intel Xeon E5-2630 v4	Dual 10-core 2.20 GHz Intel Xeon E5-2630 v4	Dual 10-core 2.20 GHz Intel Xeon Silver 4114	Dual 14-core 1.70 GHz Intel Xeon E5-2650L v4	Dual 8-core 2.40 GHz Intel Xeon E5-2630 v3	Dual 10-core 2.20 GHz Intel Xeon E5-2630 v4

Table 4. Table 4: Data describing the total time, total time spent in multigrid V-cycles, total number of multigrid V-cycles, average time spent per multigrid V-cycle, and the processor details for each simulation. This data applies to the randomly initialized simulations with Lees–Edwards boundary conditions. The number of required multigrid V-cycles decreases as the background χ 𝜒 \chi field increases, likely due to more homogeneous dynamics. Each simulation uses 32 32 32 processes.

	${\mu }_{\chi }=450\text{ K}$	${\mu }_{\chi }=500\text{ K}$	${\mu }_{\chi }=525\text{ K}$	${\mu }_{\chi }=550\text{ K}$	${\mu }_{\chi }=575\text{ K}$	${\mu }_{\chi }=600\text{ K}$
Total time (hours)	91.4188	93.9238	81.7296	94.3549	72.8395	68.4994
V-cycle time (hours)	63.5750	64.8320	51.2489	62.8451	44.4758	39.9291
# of V-cycles	34915	30467	26658	24697	22256	20495
Time/V-cycle (seconds)	6.5551	7.6606	6.9209	9.1607	7.1941	7.0137
Processor details	Dual 10-core 2.20 GHz Intel Xeon Silver 4114 v4	Dual 14-core 1.70 GHz Intel Xeon E5-2650L v4	Dual 8-core 2.40 GHz Intel Xeon E5-2630 v3	Dual 10-core 2.20 GHz Intel Xeon E5-2630 v4	Dual 14-core 1.70 GHz Intel Xeon E5-2650L v4	Dual 10-core 2.20 GHz Intel Xeon E5-2630 v4

Table 5. Table 5: Data describing the total time, total time spent in multigrid V-cycles, total number of multigrid V-cycles, average time spent per multigrid V-cycle, and processor details for each simulation. This data applies to the randomly initialized simulations with periodic boundary conditions and pure shear deformation. The number of required multigrid V-cycles decreases as the background χ 𝜒 \chi field increases, likely due to more homogeneous dynamics. Each simulation uses 32 32 32 processes.

	${\mu }_{\chi }=450\text{ K}$	${\mu }_{\chi }=500\text{ K}$	${\mu }_{\chi }=525\text{ K}$	${\mu }_{\chi }=550\text{ K}$	${\mu }_{\chi }=575\text{ K}$	${\mu }_{\chi }=600\text{ K}$
Total time (hours)	232.7589	171.3655	159.3944	129.3347	164.1213	206.7765
V-cycle time (hours)	176.7667	125.4595	115.4887	81.0119	114.6882	143.3524
# of V-cycles	66214	63871	62263	59750	57772	55149
Time/V-cycle (seconds)	9.6107	7.0734	6.6774	4.8811	7.1467	9.3577
Processor	Dual 10-core 2.20 GHz Intel Xeon E5-2630 v4	Dual 10-core 2.20 GHz Intel Xeon E5-2630 v4	Dual 10-core 2.20 GHz Intel Xeon Silver 4114 v4	Dual 16-core 2.10 GHz Intel Xeon E5-2683 v4	Dual 14-core 1.70 GHz Intel Xeon E5-2650L v4	Dual 14-core 1.70 GHz Intel Xeon E5-2650L v4

Equations

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

$$\frac{\mathcal{D}\mathbf{\sigma }}{\mathcal{D}t}=\frac{d\mathbf{\sigma }}{dt}-{\mathbf{L}}^{\mathsf{T}}\mathbf{\sigma }-\mathbf{\sigma }\mathbf{L}+\mathrm{tr}(\mathbf{L})\mathbf{\sigma },$$

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

$$\nabla \cdot \mathbf{\sigma }=0,$$

$$\frac{{\mathbf{\sigma }}^{\ast }-{\mathbf{\sigma }}^n}{\Delta t}=({\mathbf{L}}^n{)}^{\mathsf{T}}{\mathbf{\sigma }}^n+{\mathbf{\sigma }}^n{\mathbf{L}}^n-\mathrm{tr}({\mathbf{L}}^n){\mathbf{\sigma }}^n-({\mathbf{v}}^n\cdot \nabla ){\mathbf{\sigma }}^n-\mathbf{C}:({\mathbf{D}}^{\text{pl}}{)}^n.$$

$$\frac{{\mathbf{\sigma }}^{n+1}-{\mathbf{\sigma }}^{\ast }}{\Delta t}=\mathbf{C}:{\mathbf{D}}^{n+1}.$$

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

$${\tau }_0{D}^{\text{pl}}=e^{-e_z/k_B\chi }{e}^{-\Delta /k_{B}T}\cosh \left(\frac{\Omega {ϵ}_0\bar{s}}{k_{B}T}\right)\left(1-\frac{s_Y}{\bar{s}}\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).$$

$$\mathbf{x}=\mathbf{TX},$$

$$u_i=x_i-{\mathcal{R}}_i.$$

$$\mathbf{u}=\mathbf{TX}-\mathcal{R}.$$

$$\mathbf{v}=\frac{\partial \mathbf{T}}{\partial t}\mathbf{X}+\mathbf{TV}.$$

$$\mathbf{V}={\mathbf{T}}^{-1}\left(\mathbf{v}-\frac{\partial \mathbf{T}}{\partial t}\mathbf{X}\right).$$

$${\nabla }_{\mathbf{x}}={\mathbf{T}}^{-\mathsf{T}}{\nabla }_{\mathbf{X}}.$$

$$\frac{\partial \mathbf{V}}{\partial t}=-\left(\mathbf{V}\cdot {\nabla }_{\mathbf{X}}\right)\mathbf{V}+\frac{\partial {\mathbf{T}}^{-1}}{\partial t}\mathbf{TV}+{\mathbf{T}}^{-1}\left({\mathbf{T}}^{-\mathsf{T}}{\nabla }_{\mathbf{X}}\cdot \left(\mathbf{T}\mathbf{\Sigma }{\mathbf{T}}^{\mathsf{T}}\right)-\frac{{\partial }^2\mathbf{T}}{\partial t^2}\mathbf{X}-\frac{\partial \mathbf{T}}{\partial t}\mathbf{V}\right).$$

$$\mathbf{\Sigma }={\mathbf{T}}^{-1}\mathbf{\sigma }{\mathbf{T}}^{-\mathsf{T}}.$$

$$\frac{\partial \mathbf{\Sigma }}{\partial t}$$

$$\mathbf{L}={\mathbf{T}}^{-\mathsf{T}}\frac{\partial {\mathbf{T}}^{\mathsf{T}}}{\partial t}+{\mathbf{T}}^{-\mathsf{T}}{\nabla }_{\mathbf{X}}\mathbf{V}{\mathbf{T}}^{\mathsf{T}}.$$

$${\mathbf{T}}^{-1}{\mathbf{\sigma }}_0{\mathbf{T}}^{-\mathsf{T}}=\mathbf{\Sigma }-\frac{1}{3}\left({\mathbf{T}}^{-1}\mathrm{tr}\left(\mathbf{T}\mathbf{\Sigma }{\mathbf{T}}^{\mathsf{T}}\right){\mathbf{T}}^{-\mathsf{T}}\right)\mathbf{I}.$$

$${\mathbf{\sigma }}_0=\mathbf{T}{\mathbf{\Sigma }}_0{\mathbf{T}}^{\mathsf{T}}+\frac{1}{3}\left(\mathbf{T}\mathrm{tr}\left(\mathbf{\Sigma }\right){\mathbf{T}}^{\mathsf{T}}-\mathrm{tr}\left(\mathbf{T}\mathbf{\Sigma }{\mathbf{T}}^{\mathsf{T}}\right)\right)\mathbf{I}.$$

$$c_0\frac{\partial \chi }{\partial t}=-c_0\left(\mathbf{V}\cdot {\nabla }_{\mathbf{X}}\right)\chi +\frac{\left({\mathbf{D}}^{\text{pl}}:{\mathbf{\sigma }}_0\right)}{s_Y}\left({\chi }_{\infty }-\chi \right)+l^2{\mathbf{T}}^{-\mathsf{T}}{\nabla }_{\mathbf{X}}\cdot \left(D^{\text{pl}}{\mathbf{T}}^{-\mathsf{T}}{\nabla }_{\mathbf{X}}\chi \right).$$

$$\mathbf{T}=\left(\begin{array}{ccc}1& 0& U_{b}t\\ 0& 1& 0\\ 0& 0& 1\end{array}\right),$$

$$\frac{{\mathbf{\Sigma }}^{\ast }-{\mathbf{\Sigma }}^n}{\Delta t}$$

$$+{\left({\nabla }_{\mathbf{X}}\mathbf{V}\right)}^n{\mathbf{\Sigma }}^n-{\left({\mathbf{T}}^{-1}\right)}^n\mathbf{C}:{\left({\mathbf{D}}^{\text{pl}}\right)}^n{\left({\mathbf{T}}^{-\mathsf{T}}\right)}^n$$

$$\frac{{\mathbf{\Sigma }}^{n+1}-{\mathbf{\Sigma }}^{\ast }}{\Delta t}={\left({\mathbf{T}}^{-1}\right)}^n\left(\mathbf{C}:{\mathbf{D}}^{n+1}\right){\left({\mathbf{T}}^{-\mathsf{T}}\right)}^n,$$

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

$$-{\nabla }_{\mathbf{x}}\cdot {\mathbf{\sigma }}^{\ast }=-{\mathbf{T}}^n{\nabla }_{\mathbf{X}}\cdot {\mathbf{\Sigma }}^{\ast }.$$

$$\Delta t\left(\mathbf{T}^{-\mathsf{T}}\right)^{n}\nabla_{\mathbf{X}}\cdot\mathbf{C}:\bigg{(}\left(\mathbf{T}^{-\mathsf{T}}\right)^{n}\left(\nabla_{\mathbf{X}}\mathbf{V}\right)^{n+1}\left(\mathbf{T}^{\mathsf{T}}\right)^{n}\bigg{)},$$

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

$$+l^2({\mathbf{T}}^{-\mathsf{T}}{)}^n{\nabla }_{\mathbf{X}}\cdot \left((D^{\text{pl}}{)}^n({\mathbf{T}}^{-\mathsf{T}}{)}^n{\nabla }_{\mathbf{X}}{\chi }^n\right).$$

Full text

I Introduction

Molecular dynamics (MD) simulations, whereby atoms or molecules are individually simulated according to Newton’s laws [1], are widely used across the physical sciences [2, 3, 4, 5]. Open source sofware packages such as LAMMPS [6, 7] and GROMACS [8] have enabled simulations to be performed with millions of particles on modern parallel computer hardware. MD simulations provide a detailed view of the material physics and are able to capture discrete particle-level effects [9, 10]. Despite these advantages, MD simulations are computationally expensive, and it is usually only possible to simulate microscopic material samples. Furthermore, since the simulations must resolve rapid interaction timescales between particles, the applied deformation rates in MD are often orders of magnitude larger than deformation rates in laboratory tests [11, 12, 13].

Because MD simulations simulate microscopic domains, it is difficult to apply deformation via moving walls, as simulation data may be affected by finite-size effects [14, 15]. Instead, the standard approach is to apply periodic boundary conditions, but manipulate the periodic images of the primary simulation domain to achieve different applied deformations. For example, in three-dimensional Lees–Edwards boundary conditions, the periodic images have a horizontal velocity proportional to their $z$ position in order to impose simple shear [16] (Fig. 1(a)). The Kraynik–Reinelt boundary conditions [17, 18, 19, 20], plus a recent generalization by Dobson [21], use a combination of moving periodic images and domain remapping in order to simulate different extensional flows.

A complementary approach to MD is to use continuum modeling, which has the ability to simulate large system sizes and long, physically realistic timescales. However, continuum-scale theories involve a substantial theoretical hurdle, in that the transition from a particle-level theory to a continuum theory involves a coarse-graining procedure. The coarse-graining procedure defines a representative volume element (RVE) [22, 23] throughout which local deviations of material field values from their average within the RVE are neglected. The fundamental assumption of every continuum theory is that such an RVE is well-defined, and that neglecting the discrepancy between the relevant system variables and their mean within an RVE is well-justified [24, 25].

In effect, coarse-graining reduces the complex many-body system of interacting particulate constituents to a much lower degree-of-freedom system well-described by a set of nonlinear partial differential equations. This reduction in complexity is primarily responsible for the well-behaved scaling with system size in continuum simulations, in that all the classical techniques of numerical analysis become available for evolving the system over time. However, the process of coarse-graining to the continuum is difficult in general, and has primarily been successful when tailored to specific phenomena. The coarse-graining procedure introduces internal state variables that summarize the many particulate degrees of freedom, and accurate initial conditions for such internal variables can be difficult to construct. Some equilibrium systems are amenable to rigorous approaches by explicitly averaging over unwanted degrees of freedom in the system partition function [26, 27], but these approaches are intractable for many out-of-equilibrium systems.

To quantitatively explore the effect of coarse-graining MD to the continuum, it is therefore useful to perform the two types of simulation using the same geometry and conditions. However, precisely recreating the boundary conditions from MD for use in continuum simulations poses some numerical challenges. Consider the Lees–Edwards boundary conditions and suppose that the primary simulation domain is discretized on a Cartesian grid. Because the periodic images are moving, their grids will generally not align with the primary domain. This could be handled numerically via interpolation, but grid points near the boundary will incur different discretization errors. If the continuum model involves an elliptic problem, then the shifted grids will result in a complex connectivity structure in the associated linear system, which is less well-suited to some numerical linear algebra techniques.

In this work, we address this problem by developing a continuum solid mechanics simulation that permits MD boundary conditions to be recreated precisely. We use the hypo-elastoplasticity model [28] in which the deformation rate tensor $\mathbf{D}$ is decomposed additively into a sum of elastic and plastic parts [29]. There are a number of different frameworks for simulating elastoplastic materials [30], but the hypo-elastoplastic model is well-suited for problems that involve large plastic deformation. This regime is appropriate for matching to typical MD simulations, where large total strain may be applied.

Combining the additive decomposition of $\mathbf{D}$ with Newton’s second law results in a closed hyperbolic system of partial differential equations (PDEs) for the velocity $\mathbf{v}$ and stress $\bm{\sigma}$, plus coupling to evolution equations for any internal state variables. Due to the small size of MD simulations, it is usually a good approximation to say that elastic waves are fast compared to the simulation timescale, allowing for Newton’s second law to be replaced by the constraint that the stresses remain in quasi-static equilbrium, $\nabla\cdot\bm{\sigma}=\mathbf{0}$.

The resulting constrained PDE system has a mathematical correspondence to the incompressible Navier–Stokes equations, where the fluid velocity must satisfy the constraint that $\nabla\cdot\mathbf{v}=0$. For incompressible fluids a standard numerical technique is the projection method of Chorin [31, 32]. By exploiting the mathematical correspondence, a new projection method for quasi-static hypo-elastoplasticity was recently introduced [33] and extended to three dimensions [34] (Sec. II).

To match the MD boundary conditions, we introduce a coordinate transformation framework for the quasi-static hypo-elastoplastic system. It is based on an abstract linear mapping $\mathbf{T}(t)$ from a reference domain to the physical domain (Fig. 1(b)). Lees–Edwards conditions can be implemented in the continuum setting with this methodology by imposing shear through a transformation, and additionally enforcing periodic boundary conditions in all directions. Effectively, our method decouples the application of material deformation from the application of a specific boundary condition.

In addition to Lees–Edwards boundary conditions, the transformation framework is flexible, and enables simple implementation of otherwise potentially difficult deformation, such as pure shear. Any applied deformation that can be written down as a linear transformation of a reference domain can be implemented just by implementing the matrix and its time derivatives. We show that the projection method for hypo-elastoplasticity can be generalized to simulate this case by working with transformed velocities and stresses in the reference domain. The projection step in the method requires solving an elliptic problem for the velocity, and the resulting linear system has a simple mathematical structure that is well-suited for solution via numerical linear algebra techniques such as the geometric multigrid method [35, 36].

The new method is capable of simulating a wide range of elastoplastic materials, but here we consider the example of a bulk metallic glass (BMG), a new type of alloy where the atoms have a random and amorphous arrangement, in constrast to most metals [37]. BMGs have attracted considerable research interest during the past two decades. They have many favorable properties, such as high strength and wear resistance, that make them attractive for a variety of applications [38]. However, the amorphous arrangement of atoms makes the study of dynamic mechanical phenomena—such as deformation and failure—in these materials exceptionally challenging [39].

To date, a general theory of the microscopic origins of plastic deformation in amorphous solids has remained elusive. However, several prominent theories capable of making accurate qualitative and quantitative predictions have been developed, such as free-volume based theories [40, 41, 42, 43] and the shear transformation zone (STZ) theory [44, 45, 46, 47, 48]. Ultimately, free-volume theories and the STZ theory are flow-defect theories that attempt to connect microscopic rearrangements of groups of atoms with macroscopic plastic deformation, in rough analogy to the dislocation-mediated theory of plasticity in crystalline materials [49].

We employ an elastoplastic model of a BMG based on the STZ theory. A key feature of the model is the effective temperature (Sec. II.3), which characterizes the amorphous particle structure via a continuum field [50, 51, 52, 53]. The effective temperature can be measured indirectly [54], but there is currently no complete method to connect it to the microscopic particle configuration. This was recently explored by Hinkle et al. [55], who directly compared continuum and MD simulations, and examined how measurable features of MD such as the coarse-grained atomic potential energy are connected to the effective temperature. A key limitation of this study is that the MD simulations use Lees–Edwards conditions, whereas the deformation was imposed in the continuum simulation using moving parallel plates, meaning that the two could not be exactly compared. The numerical techniques that we develop here remove this limitation.

The STZ theory has proven useful for examining the failure properties of BMGs. The elastoplastic model that we employ has been used to explain the large experimental variations in notched fracture toughness of BMGs [56]. This was subsequently extended to make predictions about BMG fracture toughness for a range of parameters [57]. Recent experimental work suggests that these predictions are broadly correct [58]. BMGs also exhibit shear bands, a strain-softening instability characterized by the localization of shear strains along a thin band [59], which can be the precursor to failure [60, 61, 62]. In our simulations, we examine how shear bands nucleate as a function of the initial inhomogeneities in the effective temperature field.

The paper is organized as follows. In Sec. II, we describe the equations of quasi-static hypo-elastoplasticity and provide an introduction to the physics of the STZ theory of amorphous plasticity. In Sec. III, we introduce the coordinate transformation methodology and develop the transformed projection method. In Sec. IV we provide numerical experiments demonstrating convergence of the solution of the transformed method to the original quasi-static method in physically equivalent situations as the grid spacing is decreased. In Sec. V, we study shear banding in a bulk metallic glass subject to simple shear, Lees–Edwards, and pure shear boundary conditions. We highlight differences in results between Lees–Edwards and simple shear boundary conditions and examine how the shear band formation depends on the initial effective temperature.

II Mathematical preliminaries

II.1 Quasi-static hypo-elastoplasticity

We consider an elastoplastic material with Cauchy stress tensor $\bm{\sigma}(\mathbf{x},t)$ and velocity field $\mathbf{v}(\mathbf{x},t)$. We denote by $\mathbf{L}=\nabla\mathbf{v}$ the velocity gradient tensor and $\mathbf{D}=\frac{1}{2}\left(\mathbf{L}+\mathbf{L}^{\mathsf{T}}\right)$ the rate of deformation tensor. We adopt the framework of hypo-elastoplasticity, which assumes 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}}$. Writing linear elasticity in rate form yields

[TABLE]

where $\mathbf{C}$ is the stiffness tensor. For simplicity, the material is taken to be isotropic and homogeneous, so that $C_{ijkl}=\lambda\delta_{ij}\delta_{kl}+\mu\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right)$ where $\lambda$ is Lamé’s first parameter and $\mu$ is the shear modulus. The time derivative in Eq. 1 is the Truesdell derivative 111This expression is typically presented with $\mathbf{L}$ transposed with respect to the definition here. We adopt the convention that $\left(\nabla\mathbf{f}\right)_{ij}=\partial_{i}f_{j}$ for a vector field , i.e., partial derivatives go row-wise in gradients of vector fields. Typically, the symbol $\mathbf{L}=\nabla\mathbf{v}$ is used to denote the Jacobian or Fréchet derivative of $\mathbf{v}$, which formally is the transpose of the gradient [73]. The transformation formalism developed in this work involves both Jacobians and vector field gradients, and for physically consistent answers it is necessary to make this distinction.,

[TABLE]

with $\frac{d}{dt}=\frac{\partial}{\partial t}+\mathbf{v}\cdot\nabla$ denoting the advective derivative. The velocity field satisfies a continuum version of Newton’s second law,

[TABLE]

with $\rho$ the material density. Taken together, Eqs. 1 & 3 form a closed hyperbolic system that could form the basis of a numerical method. However, an explicit numerical method used to solve this system will resolve elastic waves. Stable resolution of elastic waves places a limit on the simulation timestep according to the well-known Courant–Friedrichs–Lewy (CFL) condition [64]. The CFL condition requires $\Delta t\leq\frac{h}{c_{e}}$ where $c_{e}$ is a typical elastic wave speed and $h$ is the grid spacing.

In metals and other materials of interest, the elastic wave speed $c_{e}$ can be large, and the grid spacing $h$ needed to resolve fine-scale features such as shear bands can be small. 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. It is often appropriate to take the long-timescale and small-velocity limit, in which the material acceleration is negligible and Eq. 3 can be replaced by the constraint

[TABLE]

which states the stresses remain in quasi-static equilibrium, and conveniently avoids the description of elastic waves. In this quasi-static limit, Eq. 1 depends on the material velocity field through $\mathbf{D}$, but the evolution equation for the velocity field has been exchanged for the constraint in Eq. 4. It is thus unclear how to solve Eq. 1 subject to the global constraint in Eq. 4.

II.2 Projection method

As noted by Rycroft et al. [33], Eqs. 1 & 3 have a close mathematical correspondence with the Navier–Stokes equations for incompressible fluid flow. The Navier–Stokes equations consist of an explicit partial differential equation for the fluid velocity along with a constraint that the velocity must be divergence-free. Much like Eqs. 1 & 4, the constraint on the velocity field is obtained from a limiting procedure applied to an explicit partial differential equation for the pressure, and the equation for the velocity still depends on the pressure after this limit has been taken.

In this setting, a well-established numerical technique is the projection method of Chorin [31, 32]. In Chorin’s projection method, the update for the velocity field is split into two steps. In the first step, an intermediate velocity field is computed which does not obey the divergence-free constraint. In the second step, a linear system is solved for the pressure field which simultaneously projects the intermediate velocity field onto the manifold of divergence-free solutions.

By using the correspondence between quasi-static hypo-elastoplaticity and incompressible fluid flow, Rycroft et al. [33] developed a new projection method for quasi-static elastoplasticity. Consider taking a timestep of size $\Delta t$, and use superscripts of $n$ and $n+1$ to denote the simulation fields before and after the step, respectively. To begin, an intermediate stress $\bm{\sigma}^{*}$ is computed by dropping the $\mathbf{C}:\mathbf{D}$ term in Eq. 1 to obtain

[TABLE]

If the velocity $\mathbf{v}^{n+1}$ were known, and hence if the total deformation rate $\mathbf{D}^{n+1}$ could be calculated, then the final stress would be given by

[TABLE]

Taking the divergence of this equation and enforcing that $\nabla\cdot\bm{\sigma}^{n+1}=\mathbf{0}$ yields

[TABLE]

After finite-difference expansion of the definition of $\mathbf{D}^{n+1}$, Eq. 7 forms a linear system for the velocity field $\mathbf{v}^{n+1}$ with source term given by the known vector $-\nabla\cdot\bm{\sigma}^{*}$, and it can be solved via standard techniques of numerical linear algebra. After solution of Eq. 7, $\bm{\sigma}^{n+1}$ can be computed according to Eq. 6, which can be shown to orthogonally project $\bm{\sigma}^{*}$ onto the manifold of quasi-static solutions. In this manner, the two-step projection method enables solving Eq. 1 subject to the global constraint Eq. 4 despite the dependence of Eq. 1 on $\mathbf{v}$. We refer the reader to papers by Rycroft et al. [33], and Rycroft and Boffi [34] for complete details on this method.

II.3 Plasticity model

As our plasticity model for a bulk metallic glass, we use an athermal form of the shear transformation zone (STZ) theory of amorphous plasticity suitable for studying glassy materials below the glass transition temperature [65, 47]. The STZ theory postulates that ephemeral and localized fluctuations of the configurational bath—STZs—occur sporadically throughout an otherwise elastic material. The STZs may be conceptualized as clusters of atoms susceptible to shear-induced configurational rearrangements when local stresses surpass the material yield stress $s_{Y}$. Each such rearrangement leads to a small increment of plastic strain, and many such rearrangements conspire to bring about macroscopic plastic deformation.

In the athermal theory used here, thermal fluctuations of the atomic configurations are neglected, and molecular rearrangements are assumed to be driven entirely by external mechanical forces. Thermal theories introduce an additional coupling between the configurational subsystem governing the rearrangements that occur at STZs, and a kinetic/vibrational subsystem governing the thermal vibrations of atoms in their cage of nearest neighbors [66]. Such thermal theories, with an additional field tracking the thermodynamic temperature that evolves according to a diffusion equation, could in principle be incorporated into our framework.

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 sits 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$ [53, 50, 51, 52].

$\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 it can be obtained as the derivative of a configurational energy with respect to a configurational entropy [39]. $\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.

We define the deviatoric stress tensor $\bm{\sigma}_{0}=\bm{\sigma}-\frac{1}{3}\operatorname{tr}(\bm{\sigma})\mathbf{I}$. The total rate of plastic deformation tensor is proportional to the deviatoric stress, $\mathbf{D}^{\text{pl}}=D^{\text{pl}}\frac{\bm{\sigma}_{0}}{\bar{s}}$, where $\bar{s}^{2}=\frac{1}{2}\bm{\sigma}_{0,ij}\bm{\sigma}_{0,ij}$ is a local scalar measure of the total deviatoric stress. The STZ theory provides the magnitude of the plastic rate of deformation as

[TABLE]

$\tau_{0}$ is a molecular vibration timescale, $e_{z}$ is a typical STZ formation energy, $k_{B}$ is the Boltzmann constant, $T$ is the thermodynamic temperature, $\Delta$ is a typical energetic barrier for a transition, $\Omega$ is a typical STZ volume, and $\epsilon_{0}$ is a typical local strain. The effective temperature satisfies a heat equation [65, 67, 68, 69, 53]

[TABLE]

The interdependence of Eqs. 8 & 9 enables the development of shear bands through a positive feedback mechanism, as increasing one of $\chi$ or $D^{\text{pl}}$ also leads to an increase in the other [67, 69].

III Coordinate transformation framework

Let $\mathbf{T}(t)$ denote a time-varying mapping from a reference domain $\mathbf{X}$ to the physical domain of interest $\mathbf{x}$ such that

[TABLE]

as shown in Fig. 1(b). Here, $\mathbf{X}\in[a_{X},b_{X}]\times[a_{Y},b_{Y}]\times[a_{Z},b_{Z}]$. Capital letters denote quantities in the reference frame and lower case letters denote quantities in the physical frame. $\nabla_{\mathbf{X}}$ and $\nabla_{\mathbf{x}}$ denote spatial differentiation in the reference and physical frame, respectively. We emphasize that $\mathbf{X}$ exists in a fixed frame on which the quasi-static hypo-elastoplastic equations will be solved, and not in the Lagrangian frame of coordinates. To clarify this point, let $\mathcal{R}=(\mathcal{X},\mathcal{Y},\mathcal{Z})$ denote a set of fixed Lagrangian coordinates. For an Eulerian frame $(x,y,z)$, we define the Eulerian displacements,

[TABLE]

We then define the Eulerian velocities $v_{i}=\frac{\partial u_{i}}{\partial t}|_{\mathcal{R}}$. The same procedure can be performed in the reference frame. We first define the physical displacements,

[TABLE]

Taking a time derivative of both sides of Eq. 12 at fixed Lagrangian coordinates $\mathcal{R}$, we arrive at an expression for the physical velocity,

[TABLE]

Above, we have identified the transformed velocity $\mathbf{V}=\frac{\partial X}{\partial t}|_{\mathcal{R}}$. Equation 13 can be used to compute the physical velocity $\mathbf{v}$ from the transformed velocity $\mathbf{V}$, if $\mathbf{V}$ is known. By inversion, it can also be used as a definition of the transformed velocity,

[TABLE]

Using the chain rule, spatial derivatives are transformed as

[TABLE]

Taking an advective time derivative of Eq. 14, using Eq. 3 for $\dot{\mathbf{v}}$, and transforming physical spatial derivatives to transformed spatial derivatives, the transformed velocity evolves according to the transformed generalization of Newton’s second law,

[TABLE]

In Eq. 16, we have rewritten the advective derivative equivalently in the reference frame, $\frac{\partial}{\partial t}+\mathbf{v}\cdot\nabla_{\mathbf{x}}=\frac{\partial}{\partial t}+\mathbf{V}\cdot\nabla_{\mathbf{X}}$. The proof of this fact for an arbitrary transformation $\mathbf{T}(t)$ is shown in Appendix A. In Eq. 16, we have also defined the transformed stress tensor via the contravariant pullback,

[TABLE]

To derive an evolution equation for $\bm{\Sigma}$, we now use the linear elastic relation in Eq. 1. Taking an advective time derivative of the relation $\bm{\sigma}=\mathbf{T}\bm{\Sigma}\mathbf{T}^{\mathsf{T}}$ and inverting, the transformed stress then obeys the transformed generalization of the linear elastic constitutive law. After expansion of the Truesdell rate,

[TABLE]

In Eq. 18, $\mathbf{D}=\frac{1}{2}\left(\mathbf{L}+\mathbf{L}^{\mathsf{T}}\right)$ refers to the physical quantity. $\mathbf{L}$ can be computed in terms of the transformed variables as

[TABLE]

$\mathbf{D}^{\text{pl}}=D^{\text{pl}}\frac{\bm{\sigma}_{0}}{\bar{s}}$ appears in Eq. 18, and its form depends on the plasticity model through the constant $D^{\text{pl}}$. As reviewed in Sec. II.1, the STZ theory provides an expression given by Eq. 8. $\mathbf{D}^{\text{pl}}$ is defined and must be computed in terms of the physical deviatoric stress $\bm{\sigma}_{0}$. In line with the definition of $\bm{\Sigma}$, we can apply the contravariant pullback to $\bm{\sigma}_{0}$ and write

[TABLE]

Using the natural definition $\bm{\Sigma}_{0}=\bm{\Sigma}-\frac{1}{3}\operatorname{tr}\left(\bm{\Sigma}\right)\mathbf{I}$ and solving for $\bm{\sigma}_{0}$, we can rewrite Eq. 20 as

[TABLE]

Equation 21 enables the computation of $\bm{\sigma}_{0}$ entirely in terms of transformed quantities. We compute $\bar{s}$ by first computing the entire tensor $\bm{\sigma}_{0}$ using Eq. 21 and then compute its Frobenius norm.

The equation for the effective temperature must also be transformed, though we do not define a transformed effective temperature. This can be accomplished by transforming the derivatives,

[TABLE]

For brevity, $\mathbf{D}^{\text{pl}}$, $\bm{\sigma}_{0}$ and $D^{\text{pl}}$ refer to the physical quantities in Eq. 22 and must be computed in terms of the transformed variables in an implementation. Transformation of the diffusive term ensures that diffusion occurs in the physical frame despite being implemented in the reference frame.

Equation 18 demonstrates that our transformation methodology leaves the Truesdell rate invariant and only affects the deformation rate term $\mathbf{C}:\left(\mathbf{D}-\mathbf{D}^{\text{pl}}\right)$. This highlights a benefit of using the Truesdell rate, as opposed to using alternative rates (e.g. Green–Naghdi or Jaumann) that employ physical approximations to achieve a simpler form. For example, the Jaumann stress rate is based upon the approximation that the effect of material stretch is much smaller than the effect of rotation, so that the Jaumann formula only involves the material spin rather than the full velocity gradient tensor. If the Jaumann rate is used in the physical frame it will not perfectly transform into the Jaumann rate in the reference frame, as neglecting stretch in the physical frame is not equivalent to neglecting stretch in the reference frame.

The transformed system of equations has connections to the principle of material frame-indifference [70, 71] which states that “the constitutive laws governing internal interactions between the parts of the system should not depend on whatever external frame of reference is used to describe them” [72]. Mathematically this is done by considering a transformation of the form $\mathbf{x}=\mathbf{R}(t)(\mathbf{X}-\mathbf{X}_{0}(t))$ where $\mathbf{X}_{0}(t)$ is a time-dependent vector and $\mathbf{R}(t)$ is a time-dependent rotation [73]. If we restrict our transformation in Eq. 10 to the case when $\mathbf{T}(t)$ is a rotation, then Eq. 18 is identical to Eq. 2, but in terms of transformed variables. The first four terms of Eq. 18 are always identical, and the final term involving $\mathbf{C}:(\mathbf{D}-\mathbf{D}^{\text{pl}})$ is also identical in this case, since the rotation matrices cancel because $\mathbf{C}$ is isotropic. Hence our coordinate transformation is consistent with material frame-indifference.

It is worth considering how the transformed system of equations differs from the original system. A particular case of interest is simple shear, given the immediate application to implementation of Lees–Edwards boundary conditions. This physical situation is described by the transformation

[TABLE]

with $U_{b}$ a boundary shear velocity. Restriction to a two-dimensional plane-strain formulation reveals that the components of Eqs. 16 & 18 retain their original form with untransformed quantities replaced by transformed quantities, in addition to several new terms proportional to powers of $U_{b}t$.

III.1 Transformed projection method

We now formulate the projection method of Sec. II.2 in the reference frame. This method enables solving for $\mathbf{V}$ and $\bm{\Sigma}$ subject to the constraint in Eq. 4. In the first step (analogous to Eq. 5), the $\mathbf{C}:\mathbf{D}$ term in Eq. 18 is neglected to compute the intermediate transformed stress $\bm{\Sigma}^{*}$,

[TABLE]

If the transformed velocity at the next timestep $\mathbf{V}^{n+1}$ were known, we could compute $\mathbf{L}^{n+1}$ via Eq. 19, compute $\mathbf{D}^{n+1}$, and complete the transformed Euler step via

[TABLE]

which is analogous to Eq. 6. To compute this correction, we need to use the physical constraint Eq. 4. Enforcing that $\nabla_{\mathbf{x}}\cdot\bm{\sigma}^{n+1}=\mathbf{0}$ leads to a linear system for $\mathbf{v}$ in the physical domain given by

[TABLE]

Because $\mathbf{T}^{-1}\bm{\sigma}^{*}\mathbf{T}^{-\mathsf{T}}=\bm{\Sigma}^{*}$ and $\nabla_{\mathbf{x}}=\mathbf{T}^{-\mathsf{T}}\nabla_{\mathbf{X}}$, the right-hand side of Eq. 26 transforms according to

[TABLE]

Similarly, the left-hand side of Eq. 26 becomes

[TABLE]

where we have omitted $\mathbf{X}$-independent terms as they will be eliminated by $\nabla_{\mathbf{X}}$. Equations 27 & 28 form a complicated linear system for the transformed velocity $\mathbf{V}^{n+1}$. The appearance of the transformation $\mathbf{T}$ in front of the gradient operator $\nabla_{\mathbf{X}}$ ensures that all mixed spatial derivatives of all components of the velocity appear in each row of Eq. 28. Equation 28 is more complex than the linear system in the original quasi-static projection method, and it is dependent on the specific form of $\mathbf{T}$. The components of Eq. 28 in the specific cases of simple shear and pure shear are shown in App. B and App. C respectively.

The update for the effective temperature is handled through an explicit forward Euler step,

[TABLE]

III.2 Numerical discretization, parallelization, and multigrid solver

The explicit update for the transformed stress Eq. 24 depends on transformed spatial derivatives of the transformed velocity through $\mathbf{L}$. Similarly, the source term in the linear system for the transformed velocity Eq. 27 depends on transformed spatial derivatives of the transformed stress. We exploit this structure through a staggered grid arrangement in the reference frame with uniform spacing $\Delta x=\Delta y=\Delta z=h$. The stress tensor $\bm{\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. Further discussion of the staggered grid arrangement can be found in [34].

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 30 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 effective temperature update in Eq. 29 is computed by expanded the divergence term,

[TABLE]

Equation 31 is computed numerically by assembling the gradient vectors $\nabla_{\mathbf{X}}\chi$ and $\nabla_{\mathbf{X}}D^{\text{pl}}$ at cell centers using the standard centered difference formula,

[TABLE]

with analogous expressions for the other directions. We also must assemble the Hessian matrix $\nabla_{\mathbf{X}}\nabla_{\mathbf{X}}\chi$ using the second derivative stencils

[TABLE]

Analogous expressions for other second derivatives are obtained through Eqs. 33 & 34 by suitable replacements. The matrix $\mathbf{T}^{-1}\mathbf{T}^{-\mathsf{T}}$ is computed from its definition.

The advective derivative in Eq. 24 must be upwinded for stability; we use the second-order essentially non-oscillatory (ENO) scheme [74]. 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. 33, the ENO derivative is defined in the $X$ direction as

[TABLE]

Above, $U_{i,j,k}$ is the $X$ component of the transformed velocity at grid point $(i,j,k)$. Equation 35 uses the curvature of $f$ to switch between an upwinded three-point derivative and a centered difference. Versions of Eq. 35 in the $Y$ and $Z$ coordinates are obtained analogously.

Despite its complexity, after spatial discretization of Eq. 28, the linear system is of the form $\mathbf{A}\mathbf{y}=\mathbf{b}$, and can be solved via standard techniques of numerical linear algebra. $\mathbf{b}$ is given in block form by the source term in Eq. 27, $\mathbf{b}_{i}=-\mathbf{T}\nabla_{\mathbf{X}}\cdot\bm{\Sigma}^{*}(\mathbf{X}_{i})$, where the index $i$ runs over all grid points. $\mathbf{y}$ is also given in block form, so that $\mathbf{y}$ contains the stacked values of $\mathbf{V}$ across all grid points. The matrix $\mathbf{A}$ is sparse, and its degree of sparsity depends on the specific discretization scheme used. In the staggered centered difference scheme described here, grid point $(i,j,k)$ is only coupled to the 27 grid points in the surrounding $3\times 3\times 3$ cube.

$\mathbf{A}$ is thus most effectively reconstructed using submatrices $\mathbf{A}^{(i,j,k)}_{(k,l,m)}$, which give the coefficients of velocity values $\mathbf{V}_{(k,l,m)}$ appearing in the linear equation for $\mathbf{V}_{(i,j,k)}$. Each matrix $\mathbf{A}^{(i,j,k)}_{(k,l,m)}$ is symmetric. With this construction, we solve Eq. 28 using a custom MPI-based parallel geometric multigrid solver; for further details of the solver, and how it interfaces with the explicit updates, the reader is referred to the non-transformed algorithm description [34]. The explicit steps for $\chi$ and $\bm{\Sigma}$ in Eqs. 24 & 29 are also parallelized using MPI and domain decomposition, with further details in the non-transformed work [34].

A highlight of the transformation methodology is its flexibility and simplicity. Implementation of new boundary conditions, so-long as they can be specified in terms of a transformation $\mathbf{T}(t)$, is only as difficult as writing the transformation down. The matrices $\mathbf{A}^{(i,j,k)}_{(k,l,m)}$ do, however, depend on the form of $\mathbf{T}(t)$, and thus they need to be derived on a transformation-by-transformation basis. Furthermore, through their dependence on $\mathbf{T}(t)$, these submatrices are time-dependent and thus need to be recomputed at each timestep.

For an arbitrary $3\times 3$ transformation with nine matrix elements, the analytical computation and hand-implementation of the $\mathbf{A}^{(i,j,k)}_{(k,l,m)}$ matrices is error-prone. To remedy this, we developed a metaprogramming scheme to auto-generate the relevant code. We used Mathematica to analytically calculate Eq. 28 in terms of arbitrary matrix elements $T_{ij}(t)$, and subsequently to replace derivatives by their finite difference equivalents. Collecting coefficients accordingly in the resulting equation gives $191$ non-zero coefficients comprising the $27$ submatrices $\mathbf{A}^{(i,j,k)}_{(k,l,m)}$. We used Python to write a skeleton file that contained function primitives for $191$ C++ functions to compute each of these coefficients individually. We then used the Mathematica function splice to fill in valid C++ code that calculates the resulting coefficients in each of these functions. Finally, we again used Python to write C++ code that calls the auto-generated C++ functions to populate the submatrices. The metaprogramming scheme only needs to be run once to generate the needed code, and it does not take any meaningful amount of time to run. In this way, the multigrid system is automatically generated at each timestep, and new simulation conditions can be immediately constructed by providing the matrix $\mathbf{T}(t)$ as a $3\times 3$ matrix class implemented in C++.

IV Numerical convergence tests

In this section, we demonstrate convergence of the transformed projection method to the non-transformed method in physically equivalent situations. In all simulations, a periodic domain in $X$ and $Y$ is considered, $-L\leq X<L$, $-L\leq Y<L$ with $L=1\text{~{}cm}$. We consider both periodic and non-periodic boundary conditions in $Z$, corresponding to domains $Z\in[-\gamma L,\gamma L)$ and $Z\in[-\gamma L,\gamma L]$, respectively. $\gamma=\frac{1}{2}$ in all simulations. We measure time in terms of the natural unit $t_{s}=L/c_{s}$ with $c_{s}=\sqrt{\mu/\rho}$ the material shear wave speed. Boundary conditions in the non-periodic case are given by

[TABLE]

Elasticity and plasticity parameters are provided in Table 1, and for these parameters, $t_{s}=4.05\text{~{}\textmu{}s}$. All simulations in this section are run with $32$ processes on an Ubuntu Linux computer with dual 14-core 1.70 GHz Intel Xeon E5-2650L v4 processors.

The global three-dimensional grid has spacing $h$ in each direction. The cell-cornered grid points are indexed according to $i=0,\ldots,Q-1$, $j=0,\ldots,M-1$ in the $X$ and $Y$ directions. In the $Z$ direction, the grid points are indexed according to $k=0,\ldots,N$ and $k=0,\ldots,N-1$ for non-periodic and periodic boundary conditions, respectively. The cell-centered grid points run according to $i=\frac{1}{2},\frac{3}{2},\ldots Q-\frac{1}{2}$, $j=\frac{1}{2},\frac{3}{2},\ldots M-\frac{1}{2}$, and $k=\frac{1}{2},\frac{3}{2},\ldots M-\frac{1}{2}$. As described in Sec. III.2, $\bm{\Sigma}$ and $\chi$ are stored at cell centers while $\mathbf{V}$ is stored at cell corners. The additional grid points $(i,j,k=N)$ in the $Z$ direction in the non-periodic case are ghost points used for enforcing the Dirichlet boundary conditions $\mathbf{V}=\mathbf{0}$.

The cell-centered grid points on the top boundary $(i,j,N+\frac{1}{2})$ contain linearly extrapolated $\bm{\Sigma}$ and $\chi$ values to ensure that $\bm{\Sigma}$ and $\chi$ remain free on the top boundary. In the periodic case, the grid points $(i,j,k=N)$ hold the velocity values $\mathbf{V}_{(i,j,0)}$, and the corresponding cell-centered grid points are used to hold the wrapped values of $\bm{\Sigma}_{(i,j,\frac{1}{2})}$ and $\chi_{(i,j,\frac{1}{2})}$. At the simulation boundaries in the $X$ and $Y$ directions, ghost points leaving the simulation domain are filled with values that wrap around, so that the ghost point corresponding to grid point $(Q,j,k)$ is filled with the real values from grid point $(0,j,k)$. Similarly, values at points $(i,M,k)$ are filled using values from $(i,0,k)$.

IV.1 Qualitative comparison between the transformed and non-transformed methods

We now demonstrate the qualitative similarity of solutions computed with the transformed and the standard quasi-static methods. In the following subsection, this comparison is made quantitatively rigorous. To visualize the results three-dimensionally, we use a custom opacity function,

[TABLE]

where $\chi_{bg}$ is a background $\chi$ value. By choice of $a$ and $\eta$, the most physically relevant features in three-dimensional visualizations of the $\chi$ field can be revealed.

To compare the transformed and non-transformed methods, a physically equivalent situation is now constructed. We employ non-periodic Dirichlet boundary conditions in the $Z$-direction and enforce $\mathbf{V}(X,Y,\pm\gamma L)=(0,0,0)$. To impose deformation, we use a shear transformation $\mathbf{T}(t)$ corresponding to

[TABLE]

Boundary conditions in the non-transformed simulation correspond to shearing between two parallel plates, $\mathbf{v}(x,y,\pm\gamma L)=(U_{b},0,0)$. An initial linear velocity gradient is imposed in the non-transformed frame, so that

[TABLE]

Equation 39 ensures equivalent initial conditions in both methodologies, and also prevents the introduction of large gradients in the deformation rate near the boundary. To create interesting dynamics, an initial condition in $\chi$ corresponding to a helix oriented perpendicular to the direction of shear is considered. This is represented as

[TABLE]

Equation 40 is written in the non-transformed simulation, but the same initial conditions are used in the transformed simulation with the substitution $\mathbf{x}\rightarrow\mathbf{X}$. The configuration is visualized in the reference frame in Fig. 2.

The simulations are conducted on two grids of size $256\times 256\times 128$ with a quasi-static timestep $\Delta t=31.25t_{s}$ and with a value of $U_{b}=10^{-7}\frac{L}{t_{s}}$. Snapshots at three representative time points are shown in Fig. 3. In Fig. 3 (a,b) at $t=2.88\times 10^{5}t_{s}$, shear band nucleation has not begun, and there is an increase in the $\chi$ field across the entire domain. At $t=4.02\times 10^{5}$ in Fig. 3 (c,d), shear bands have begun to nucleate along the top and bottom planes of the helices. At $t=6\times 10^{5}t_{s}$ in Fig. 3(e,f), the bands have grown sharper, stronger, and span the system. In all cases, the qualitative agreement is very good.

IV.2 Quantitative comparison between the transformed and non-transformed methods

Having demonstrated the qualitative similarity between the solutions computed by the transformed and non-transformed methods, we now present a rigorous quantitative comparison. We utilize the same simulation geometry, boundary conditions, shear transformation, and initial conditions as in Sec. IV.1. We introduce a norm over simulation fields,

[TABLE]

where the integral in Eq. 41 runs over the entire simulation domain and is numerically computed using the trapezoid rule. The appearance of $|\cdot|$ in Eq. 41 is interpreted as the two-norm for vectors, absolute value for scalars, and the Frobenius norm for matrices. With subscript NT denoting “non-transformed” and subscript T denoting “transformed”, Eq. 41 is applied to the quantities $\mathbf{v}(\mathbf{X},t)_{\text{NT}}-\mathbf{v}(\mathbf{X},t)_{\text{T}}$, $\bm{\sigma}(\mathbf{X},t)_{\text{NT}}-\bm{\sigma}(\mathbf{X},t)_{\text{T}}$ and $\chi(\mathbf{X},t)_{\text{NT}}-\chi(\mathbf{X},t)_{\text{T}}$. The physical field values are compared across the reference grid, a procedure that involves two subtleties.

In the transformed simulation, this comparison requires computing $\bm{\sigma}$ from $\bm{\Sigma}$ and $\mathbf{v}$ from $\mathbf{V}$ using Eqs. 17 and 14 respectively at all reference grid points. In the non-transformed simulation, it is necessary to compute the non-transformed field values at reference grid points. Because the reference grid maps to a sheared physical grid, these values may not be defined in the non-transformed simulation. We handle this via the following procedure. The non-transformed simulation grid point $\mathbf{x}(\mathbf{X})$ corresponding to the reference grid point $\mathbf{X}$ is first computed. If $\mathbf{x}(\mathbf{X})$ does not lie on the non-transformed grid, adjacent grid points are linearly interpolated to compute an approximate field value at $\mathbf{x}$. This procedure incurs an $\mathcal{O}(h^{2})$ error, which is the same order of accuracy as the centered differences used for spatial discretization in the two methods. As the sizes of the simulation grids are increased, the discrepancy in solutions will decrease.

To ensure that issues with temporal discretization do not affect the comparison, it is also necessary to scale the quasi-static timestep as the grid size is decreased. Because the spatial order of accuracy is $\mathcal{O}(h^{2})$, we keep the ratio $\Delta t/h^{2}$ fixed across all simulations. We perform comparisons across grids of size $N\times N\times\frac{N}{2}$ with $N=64$, $96$, $128$, $160$, $192$, and $256$. Respectively, these correspond to grid spacings $L/32$, $L/48$, $L/64$, $L/80$, $L/96$, and $L/128$. The quasi-static timestep is taken to be $\Delta t=500t_{s}$ for the coarsest simulation, leading to quasi-static timesteps $\Delta t=222.14,125,80,55.55$, and $31.25$ respectively for the finer simulations. The diffusion length scale in the effective temperature equation is taken to be zero in all simulations for the purpose of the comparison.

The results for the quantitative comparisons are shown in Fig. 4. In Fig. 4 (top left), the three $L_{2}$ norm curves are plotted together for a value of $N=256$, where each curve is normalized by a representative value in order to plot on a comparable dimensionless scale. The effective temperature norm increases rapidly early on in the simulation, but then saturates around $10^{-4}$. The $\bm{\sigma}$ norm stays around machine precision until the onset of plasticity, when it rapidly increases and then saturates around $10^{-3}$. Similarly, the $\mathbf{v}$ norm stays below $10^{-13}$ until the onset of plasticity, when it rapidly increases and then saturates around $10^{-4}$. The agreement up to machine precision prior to the onset of plasticity is expected, and validates the accuracy of the derivation of the equations in the reference frame.

In Fig. 4 (top right), the effective temperature norm curves are shown for all values of $N$. Here, there is a steady increase in the discrepancy before the onset of plasticity due to advection across the grid. After plasticity is activated around $t=1.2\times 10^{5}t_{s}$, there is a period of saturation in all curves, followed by a period of increase beginning around $t\approx 3\times 10^{5}t_{s}$, where some simulation curves cross and end at roughly equal values. As expected, the discrepancies generally decrease as the grid spacing is decreased.

In Fig. 4 (bottom left), the velocity norm curves are shown as a function of time for all discretization levels. In all cases, the difference between the simulation methods is on the order of machine precision until the onset of plasticity, when there is a sharp and immediate jump. The size of the jump decreases with the discretization level as expected.

In Fig. 4 (bottom right), the stress norm curves are shown. These curves display a combination of the trends in the velocity and effective temperature plots. Before the onset of plasticity, the error in all simulations is very low - on the order of machine precision. After the onset of plasticity, there is a sharp jump in all simulations, and the size of the jump decreases with higher resolution. Past around $t\approx 2\times 10^{5}t_{s}$, the curves begin to cross, all ending at roughly equivalent values.

To compare the computational efficiency of the two methods, we have reported timing statistics for all simulations in Tab. 2. Displayed are the total time, the total number of multigrid V-cycles, the total time spent in multigrid V-cycles, and the average time per V-cycle for the transformed (T) and non-transformed (NT) methods. As is clear from the table, the transformed method incurs a mild increase in computational expensive. The average time spent per V-cycle is roughly the same, but the total number of multigrid V-cycles is higher for the transformed method. This is likely due to the increased complexity of the linear system required for the stress projection in the transformed method when compared to the non-transformed method.

V Numerical examples

V.1 Simple shear and the effect of Lees–Edwards boundary conditions

As a first example application of the transformation method, we consider connecting a continuum-scale model to typical discrete molecular dynamics simulations. A significant difference between continuum simulation and molecular dynamics is in the boundary conditions. Molecular dynamics simulations commonly employ Lees–Edwards boundary conditions, where periodic copies of the system are placed above and below with a prescribed horizontal velocity. Continuum-scale boundary conditions usually set a shear velocity on the top and bottom boundaries to achieve the same effect.

Lees–Edwards boundary conditions can be implemented in the continuum through the use of the coordinate transformation methodology presented here, by combining a shear transformation $\mathbf{T}(t)$ as in Eq. 23 with periodicity in the $Z$ direction. In the following sections, we present several numerical examples using Lees–Edwards and non-periodic boundary conditions. Particular attention is paid to differences in shear banding dynamics produced by these two choices of boundary conditions.

V.1.1 Cylindrical inclusion

We first consider an initial condition corresponding to a cylindrical defect in the material. Accordingly, the effective temperature field is initially elevated throughout a cylinder of finite length oriented along the direction of shear,

[TABLE]

The initial condition is shown in Fig. 5. The diffusion lengthscale is set to $l=\frac{3}{2}h$ and the quasi-static timestep is set to $\Delta t=200t_{s}$. The grid is of size $256\times 256\times 128$. The simulation is performed to a final value of $t=2\times 10^{6}t_{s}$. To induce shear banding, a shear transformation of the form Eq. 23 is used with a value of $U_{b}=10^{-7}L/t_{s}$, and both clamped and Lees–Edwards boundary conditions are considered. The clamped simulation takes 13.851 total hours when run with 32 processes on an Ubuntu Linux computer with dual 10-core 2.20 GHz Intel Xeon E5-2630 v4 processors. 10.452 hours are spent in multigrid V-cycles and 28293 total V-cycles are required. The Lees–Edwards simulation takes 10.082 total hours when run with 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 7.393 hours and 28293 total V-cycles are required.

Results for Lees–Edwards and non-periodic boundary conditions are shown in Fig. 6, on the right and left respectively. The shear banding dynamics in this case are simple, and correspond to outward nucleation of a single band from the localized cylinder. At $t=5\times 10^{5}t_{s}$ in Fig. 6 (a,b), nucleation of the shear band has begun, and there is some spreading in the $\chi$ field visible at the caps of the cylinder. By $t=1.25\times 10^{6}t_{s}$, a prominent system-spanning shear band has formed, as displayed in Fig. 6 (c,d). In Fig. 6 (e,f) at $t=2\times 10^{6}t_{s}$, the shear band continues to grow stronger and thicker. In this case, the dynamics are virtually identical for the Lees–Edwards and nonperiodic boundary conditions.

V.1.2 A randomly fluctuating effective temperature field

We now consider a randomly distributed initial condition in the effective temperature field $\chi(\mathbf{X},t=0)$. 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 extended set of grid points and ghost points. Equation 43 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. 43 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. For computational feasibility, we choose a cutoff length of $5l_{c}$, so that the Gaussian kernel is considered to be zero past this point. In the following studies, a value of $l_{c}=5h$ is used, leading to an additional $25$ ghost points padding the grid for the purpose of the convolution.

Simulations are performed for mean values $\mu_{\chi}=450\text{~{}K},500\text{~{}K},525\text{~{}K},550\text{~{}K},575\text{~{}K},600\text{~{}K}$ with a fixed value of $\sigma_{\chi}=15\text{~{}K}$ for both non-periodic and Lees–Edwards boundary conditions. The diffusion length scale is set to $l=\frac{3}{2}h$, and the quasi-static timestep is set to $\Delta t=200t_{s}$. The simulations are all conducted on a $512\times 512\times 256$ cell grid to a final value of $t=1\times 10^{6}t_{s}$. To induce shear banding, a shear transformation of the form Eq. 23 with a value of $U_{b}=10^{-7}L/t_{s}$ is imposed on the domain. Timing details for the non-periodic simulations are shown in Table 3, while timing details for the Lees–Edwards simulations are shown in Table 4.

The results for this sequence of simulations in the case of non-periodic boundary conditions are shown in Figs. 7–10. Each figure corresponds to a single snapshot in time, and the mean increases with the alphabetical labeling. The value of $\chi_{bg}$ used in the opacity function in each case is given by $\mu_{\chi}-25\text{~{}K}$. The initial conditions for the effective temperature field are shown in Fig. 7. At $t=0$, all simulations look essentially the same. The realization of the noise in each configuration is identical, and each pane is obtained from the previous by a constant shift in $\chi$.

By $t=4\times 10^{5}t_{s}$ in Fig. 8, the simulations with the two lowest values of $\mu_{\chi}$ exhibit clear shear bands with curvature in both the $X$ and $Y$ directions. The simulation with $\mu_{\chi}=450\text{~{}K}$. in Fig. 8(a) displays two shear bands that cross each other diagonally near $\frac{X}{L}=0.5$. The simulation with $\mu_{\chi}=500\text{~{}K}$ in Fig. 8(b) displays only one of these bands, though the second has begun to nucleate. This single band is also apparent in Fig. 8(c), but it is significantly weaker. A third nascent band near $\frac{Z}{L}=0$ may also be observed.

More details are clear at $t=6\times 10^{5}t_{s}$ in Fig. 9. Figure 9(a) is similar to Fig. 8(a), whereas Fig. 9(b) shows further development of the shear bands in Fig. 8(b). Figures 9(c) and (d) show the development of several flat and thin shear bands centered around $\frac{Z}{L}=0$.

Figure 10 ($t=10^{6}t_{s}$) displays clear shear banding across all values of $\mu_{\chi}$, and makes clear the dependence of shear banding structure on $\mu_{\chi}$. There is one primary band in Fig. 10(a), with a split near around $\frac{X}{L}\approx-0.5$ not present in Fig. 10(b). Figure 10(b) also displays an additional thin band near $\frac{Z}{L}\approx 0.0$ that has not formed in Fig. 10(a). Figure 10(c) displays several additional bands near $\frac{Z}{L}=0$ that form a complex branching pattern. Figure 10(d) resolves more fine-scale structure in the band near $\frac{Z}{L}\approx 0.25$ when compared to Figs. 10(a)-(c) as if a single band has begun to split, and has more bands near lower values of $\frac{Z}{L}$. Figs. 10(e) and (f) show several additional thin bands when compared to the previous panels, but they are earlier in their formation and less prominently displayed.

Taken together, Figs. 7–10 provide qualitative insight into how macroscopic shear banding dynamics and structure reflect the underlying effective temperature distribution. In a simulation with small mean, there are few regions susceptible shear band nucleation, most clearly displayed in the formation of only a single band in the lowest mean simulation. These nucleation points must connect to form a band, as indicated by the mild curvature seen in the bands in Figs. 10(a) and (b). As $\mu_{\chi}$ is increased, additional regions of sufficiently high $\chi$ exist for band nucleation, curvature decreases, and the number of bands increases. This first presents itself, as seen in Figs. 7(d)–10(d), as an existing band splitting into multiple. The gap in the split grows with $\mu_{\chi}$, as seen in Figs. 7(d)-10(d) and Figs. 7(e)-10(e), until it eventually breaks off into its own band. With high $\mu_{\chi}$ as in Figs. 7(e)-10(e) and Figs. 7(f)-10(f), shear bands can nucleate in many different locations without curvature. The timescale for shear band development is also more rapid in simulations with low background $\chi$ field.

The results for an identical sequence of simulations in the case of Lees–Edwards boundary conditions are displayed in Figs. 11–14. The initial conditions are displayed in Fig. 11, which differ from those in Fig. 7, as the convolution used to generate the initial distribution wraps around over the boundary in $Z$ to enforce periodicity.

By $t=4\times 10^{5}t_{s}$ in Fig. 12(a), a single vertical shear band has formed, along with an additional, weaker vertical band and a similar horizontal band. These bands are also visible in Fig. 12(b) earlier in their development. Vertical shear bands do not typically form in continuum simulations with non-periodic boundary conditions in $Z$, but are frequently seen in MD simulations [75, 76], indicating that the orientation of shear bands could be strongly related to the specific boundary conditions used.

Further progression is clear at $t=6\times 10^{5}t_{s}$ in Fig. 13. Figure 13(a) is similar to Fig. 12(a). Figure 13(b) displays significant strengthening of the early-stage bands in Fig. 12(b). Figs. 13(c) and (d) show the initiation of several shear bands.

Figure 14 shows the results for $t=10^{6}t_{s}$. Figure 14(a) displays a horizontal band perpendicular to the vertical band that exhibits significant curvature. Figure 14(b) shows a similar result, with a thinner vertical band and a thicker, flatter horizontal band. Figure 14(c) shows similar features, but also displays the development of several additional horizontal bands extending to the bottom of the simulation domain. Furthermore, the thick horizontal band in Fig. 14(b) can be seen to split and fracture. In Fig. 14(d), the vertical band has been almost entirely washed out, and a complex branching pattern of horizontal bands is seen. Figs. 14(e) and (f) are similar to Fig. 14(d), but they are earlier in their development and some of the fine-scale features are washed out due to the high background $\chi$ field. The agreement with the non-periodic simulations increases strongly as $\mu_{\chi}$ is increased.

Figure 14 clearly demonstrates the effect of increasing $\mu_{\chi}$ with periodic boundary conditions. In the simulations with lower $\mu_{\chi}$, nucleation of vertical shear bands is more likely, and curved horizontal bands develop later in the simulation than vertical bands. As $\mu_{\chi}$ is increased, the vertical bands begin to disappear. As in the non-periodic case, the curvature in the horizontal bands decreases with $\mu_{\chi}$. As $\mu_{\chi}$ is increased further, the vertical bands disappear altogether. In this regime, increasing $\mu_{\chi}$ increases the number of horizontal bands, and the qualitative agreement with the nonperiodic results is good. These results suggest that, for higher $\mu_{\chi}$, the effect of periodicity in the $Z$ direction is less significant.

V.2 Pure shear

As a second example transformation, we now consider pure shear deformation. In metallic glasses, experimental evidence indicates that pure shear is the primary failure mode under compressive stress, and several recent experiments have been conducted probing BMGs under pure shear conditions [77, 78, 79, 80]. Pure shear is particularly interesting due to the simplicity of its implementation in the transformation methodology. To simulate pure shear on a physical grid, it would be necessary to impose traction boundary conditions on the top, bottom, and sides, which poses computational difficulties. Within the transformation framework, pure shear can be implemented using the transformation

[TABLE]

$A(t)$ can be chosen as any monotonically increasing function of time. In the following studies, we choose $A(t)=e^{\xi t}$, where $\xi$ is a simulation parameter that sets the rate of extension and compression of the $x$ and $z$ axes respectively. For numerical stability, it is important to choose $\xi$ small, so that large stresses do not cause divergences in the simulation fields. In our simulations, we choose $\xi$ as a fraction of $t_{f}$, which effectively sets the strain at the end of the simulation.

V.2.1 Gaussian defects

To gain some physical intuition about shear banding dynamics with pure shear boundary conditions, we first consider an example initial condition in $\chi$ corresponding to localized defects in the material. It is expected that diagonal shear bands will nucleate outwards from the imperfections. We first define the quantities

[TABLE]

and then take the initial condition in $\chi$ to be

[TABLE]

Simulations are performed with periodic and non-periodic boundary conditions in $Z$ on grids of size $256\times 256\times 128$. The $X$ and $Y$ dimensions use periodic boundary conditions in both cases. The diffusion length scale is set to $l=3h$ and the quasi-static timestep is $\Delta t=200t_{s}$. $\xi$ in Eq. 44 is set to be $\frac{1}{4t_{f}}$ with $t_{f}=4\times 10^{5}t_{s}$ the total simulation duration, so that $A(t_{f})=e^{1/4}\approx 1.284$. Both periodic and non-periodic boundary conditions are considered. The periodic simulation is run with 32 processes on an Ubuntu Linux computer with dual 14-core 1.70 GHz Intel Xeon E5-2650L processors. The total time is 4.721 hours, the total time spent in multigrid V-cycles is 3.283, and the total number of V-cycles is 11596. The non-periodic simulation is run with 32 processes on an Ubuntu Linux computer with dual 10-core 2.20 GHz Intel Xeon E5-2630 processors. The total time is 8.832 hours, the total time spent in multigrid V-cycles is 6.506 hours, and the total number of V-cycles is 11075.

Results for periodic and non-periodic boundary conditions are shown in Figs. 15 and 16 respectively. The initial conditions are shown in Figs. 15(a) and 16(a). In both Figs. 15(b) and 16(b) at $t=5\times 10^{4}t_{s}$, some spreading in the $\chi$ field is seen near the defects. Shortly thereafter, the dynamics in the nonperiodic and periodic cases begin to differ dramatically.

At $t=8\times 10^{4}t_{s}$ in Fig. 15(c), three diagonal bands are seen connecting the defects. The bands become more pronounced at $t=10^{5}t_{s}$ in Fig 15(d). This continues into $t=2\times 10^{5}t_{s}$ in Fig. 15(e), along with the addition of diagonal bands perpendicular to the original bands. Both bands continue to grow larger and stronger by $t=4t\times 10^{5}t_{s}$ in Fig. 15(f).

The deformation dynamics with non-periodic boundary conditions are significantly different. By $t=8\times 10^{4}t_{s}$ in Fig. 16(c), diagonal bands have started to nucleate off of each defect in a direction roughly perpendicular to the first bands formed in the periodic simulation. By $t=10^{5}t_{s}$ in Fig. 16(d), this nucleation has grown more prominent, and an increase in the background $\chi$ field is seen across the simulation. At times $t=2\times 10^{5}t_{s}$ and $t=2.5\times 10^{5}t_{s}$ in Figs. 16(e) and 16(f) respectively, the qualitative structure remains the same, but the background $\chi$ field continues to increase. Unlike in the periodic case, true system-spanning shear bands do not fully form.

V.2.2 A randomly fluctuating effective temperature field

In this section, we consider the same sequence of random initializations as in Sec. V.1.2, but now subject to pure shear deformation. Simulations are performed across values of $\mu_{\chi}=450\text{~{}K}$, $500\text{~{}K}$, $525\text{~{}K}$, $550\text{~{}K}$, $575\text{~{}K}$, and $600\text{~{}K}$ with a fixed value of $\sigma_{\chi}=15\text{~{}K}$. The diffusion length scale is set to $\frac{3}{2}h$ and the quasi-static timestep is set to $\Delta t=200t_{s}$. All simulations are conducted on a $512\times 512\times 256$ cell grid. A pure shear transformation of the form Eq. 44 is used with $A(t)=e^{\xi t}$ and a value of $\xi=\frac{1}{4t_{f}}$ with $t_{f}$ = $2\times 10^{6}t_{s}$ so that $A(t_{f})=e^{1/4}\approx 1.284$. Simulations are performed with fully periodic boundary conditions in all directions; non-periodic simulations produce qualitatively similar differences as in the case of simple shear. In all figure panels, $\chi_{bg}$ is set to be $\mu_{\chi}-25\text{~{}K}$. Timing data for the simulations is reported in Table 5.

The results are shown in Figs. 17–20, with the initial condition shown in Fig. 17. All simulations undergo an increase in $\chi$ until the formation of diagonal shear bands begins. Much like the defect simulations seen in the previous section, shear bands nucleate diagonally at roughly $45\degree$ angles to the $X$–$Y$ plane. As in the simple shear simulations, distributions in $\chi$ with higher mean values have slower dynamics. The structural effect of varying $\mu_{\chi}$ is most easily seen in Fig. 20. As $\mu_{\chi}$ increases, the number of shear bands vastly increases, forming a cross-hatched pattern throughout the domain. The cross-hatching becomes more regular and more finely spaced with higher values of $\mu_{\chi}$.

VI Conclusion

In this work, we derived the equations of hypo-elastoplasticity on a fixed reference domain which can be mapped to a physically deforming material through a time-varying linear transformation $\mathbf{T}(t)$. The difference between this frame and the Lagrangian frame was shown, and the utility of this frame in implementing complex boundary conditions such as the Lees–Edwards conditions used in molecular dynamics and pure shear in a fully periodic setting was demonstrated. The quasi-static projection algorithm was derived in the reference frame and its convergence to the standard method was shown as the level of discretization increases. Several numerical examples were considered in the STZ model of amorphous plasticity. In particular, for a randomly-distributed initial condition in the effective temperature field, the dependence of shear banding dynamics on the mean of the distribution was discussed under conditions of simple shear and pure shear. Our work also highlights that the direction of shear bands (e.g. horizontal versus vertical in simple shear) can be strongly influenced by boundary conditions.

With the simple implementation of Lees–Edwards conditions afforded by the transformation method, boundary conditions can now be made equivalent in MD and continuum modeling. The development of a method to compute a precise matching between atomic configurations in molecular dynamics and effective temperature distributions in continuum simulations is a promising direction of future research. The ability to do so would place internal state variables in plasticity models (such as the effective temperature in the STZ model) on a firmer theoretical footing. In addition, hybrid computational approaches could be developed, where an MD simulation could first be used to compute an initial condition for a significantly larger scale continuum simulation. This type of approach would combine the physical accuracy of MD with the capability of continuum simulations to simulate large system sizes and long times. As an added benefit, our approach enables the study of the effect of periodic boundary conditions in general, independent of the relevance of these settings to MD.

So far, our implementations are restricted to cases where the material fills the entire computational domain, and loading is applied via planar boundary conditions, or via the coordinate transformation framework. However, the methods presented here could be generalized to materials with free boundaries, using the level set method [81, 82] to track the material boundary. Methods to do this have already been implemented in two dimensions [56, 83, 33], and the same methods could be used in principle in three dimensions. However, it is a challenging computational task, since it requires extensive modifications to the finite-difference stencils near the material boundary. In particular, since some grid points will lie outside the material, the geometric multigrid method is no longer well-suited for solving the projection step, since it relies on a regular arrangement of grid points. It may be necessary to use algebraic multigrid approaches or Krylov-based linear solvers. Nevertheless this remains a high priority for future work, since it would open up many new directions, such as studying three-dimensional cavitation [84, 11], simulating mode III fracture [85], and predicting the topography of fracture surfaces [86, 87, 88].

Appendix A Advective derivative calculation

Consider a scalar field $\phi(\mathbf{x},t)=\phi(\mathbf{T}\mathbf{X},t)$. We can compute the advective derivative of $\phi$ as follows using the chain rule,

[TABLE]

In the last line, we have used Eq. 14 and the identity $\frac{\partial\mathbf{T}^{-1}}{\partial t}=-\mathbf{T}^{-1}\frac{\partial\mathbf{T}}{\partial t}\mathbf{T}^{-1}$.

Appendix B Linear system for simple shear

Let $\mathbf{V}=\left(U,V,W\right)^{\mathsf{T}}$. For a simple shear transformation as given in Eq. 23, $\mathbf{C}:\mathbf{D}$ takes the form

[TABLE]

The above set of equations leads to the linear system (Eqs. 27 & 28) for the velocity field

[TABLE]

Discretization of the second-derivative terms in Eqs. B–55 using the finite differences in Sec. III.2 enables application of the geometric multigrid method to solve for $U$, $V$, and $W$.

Appendix C Linear system for pure shear

For a pure shear transformation of the form Eq. 44 with $A(t)=e^{\xi t}$ as in the main text, $\mathbf{C}:\mathbf{D}$ takes the form

[TABLE]

The above set of equations leads to the linear system (Eqs. 27 & 28) for the velocity field

[TABLE]

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. 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.