A multiscale shear-transformation-zone (STZ) model and simulation of plasticity in amorphous solids

Accepted Manuscript A multiscale shear-transformation-zone (STZ) model and simulation of plasticity in amorphous solids Shingo Urata, Shaofan Li PII:

S1359-6454(18)30425-7

DOI:

10.1016/j.actamat.2018.05.058

Reference:

AM 14611

To appear in:

Acta Materialia

Received Date: 26 January 2018 Revised Date:

21 April 2018

Accepted Date: 25 May 2018

Please cite this article as: S. Urata, S. Li, A Multiscale Shear-Transformation-Zone (STZ) Model and Simulation of Plasticity in Amorphous Solids, Acta Materialia (2018), doi: 10.1016/ j.actamat.2018.05.058. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

A Multiscale Shear-Transformation-Zone (STZ) Model and Simulation of Plasticity in Amorphous Solids Shingo Urata † , Shaofan Li †

‡∗

RI PT

Innovative Technology Research Center, Asahi Glass Co. Ltd. 1-1 Suehiro-cho, Tsurumi-ku, Yokohama, Kanagawa, 230-0045, Japan ‡

SC

Department of Civil and Environmental Engineering, The University of California, Berkeley, CA 94720, USA

Abstract

EP

TE

D

M AN U

In this work, we have developed a multiscale shear-transformation-zone (STZ) model to simulate inelastic deformations in amorphous materials at macroscale without using phenomenological constitutive modeling. The main novelties of the multiscale method are two: (1) the multiscale method employs an atomistic-based representative sampling cell (RS-cell) that is embedded in every quadrature points of macroscale finite elements, and it generically contains potential shear-transformation-zones. Thus the evolution of RS-cells can naturally allow molecular clusters having irreversible microstructure rearrangements at microscale in response to applied loads, and subsequently translate the effects of the irreversible rearrangements of molecular clusters in each RS-cell to inelastic responses at the continuum scale, i.e. delineating the defect-mediated amorphous plasticity, and (2) the method uses a ParrinelloRahman molecular mechanics based the Cauchy-Born rule to construct an atomisticallyinformed constitutive model at continuum level, which is able to quantitatively measure amorphous plasticity behaviors, including plastic flow stress, inelastic hysteresis loops under cyclic loading, and strain localizations at macroscale. By doing so, we have successfully simulated plastic deformation in the Lennard-Jones binary glass (LJBG) at macroscale level. Moreover, we have shown that the multiscale shear-transformation-zone can capture cyclic plasticity and the shear band formations in LJBG material.

AC C

Keywords: Amorphous materials, Molecular dynamics, Cauchy-Born rule, Multiscale simulation, Parrinello-Rahman molecular dynamics, Plasticity

1. Introduction

Structurally discorded amorphous materials are noncrystalline solids, and they include metallic glasses, glassy polymers, amorphous ceramics, amorphous thin films, some granular ∗

Corresponding Author Email address: [email protected] (Shaofan Li ‡ )

Preprint submitted to Acta Materialia

June 5, 2018

ACCEPTED MANUSCRIPT

AC C

EP

TE

D

M AN U

SC

RI PT

materials, and many biological materials, etc. In past fifty more years, many researches have been conducted to study the mechanical properties of amorphous materials, because of their broader applications ranging from construction materials, electronic materials, energy materials, to bio-materials, in particular their potential to become a class of mainstream structural materials. Many amorphous materials are quasi-brittle materials that have high elastic strength with some ductility. To understand the origin of such ductility as well as failure mechanism will be instrumental for controlling the material responses and enhance its performance. Due to their disordered and random microstructures, mechanical behaviors of amorphous materials at macroscale are difficult to predict by using the modeling approach in comparison with that of crystalline materials, because their random microstructure may lead to a complex energy landscape at microscale, and a slight change of deformation field may lead to an irreversible change of the minimum energy path in the free energy surface. In other words, atom positions at microscale can be deformation history dependent, even when the whole amorphous system is still in a macroscopically elastic state. In recent years, the field has advanced significantly, and this is partially due to the insight gained from molecular dynamics based simulations, which have inspired new theoretical developments, and one of the prominent plasticity theories in amorphous materials is the shear-transformation-zone (STZ) theory. There are mainly two versions of STZ theory: (I) a defect-mediated plasticity theory e.g. Spaepen [1], Argon and his co-workers [2, 3], and more recently Homer and Schuh [4], and (II) the non-equilibrium thermodynamics based STZ theory developed by Falk and Langer and their co-workers, e.g. [5, 6, 7, 8, 9, 10], which places its emphasis on the non-equilibrium thermodynamics origin and foundation of plasticity of disordered materials. In this work, the term of shear-transformation-zone is mainly referred to the former, in which localized clusters of aggregated molecules undergo irreversible rearrangements in response to applied shear stresses. In recent years, there have been a large body of literatures published on the molecular dynamics foundation of the defect-mediated STZ theory as well as its applications, in particular on metallic glasses, e.g. [5, 11, 12, 13, 14, 15, 16]. On the other hand, however, the STZ theory is a phenomenological defect-mediated flow theory, and it is not a microstructure-specific theory, e.g. it does not take into account the specific disordered microstructure in different amorphous materials, which may result in different inelastic behaviors. Moreover, defects in disordered noncrystalline solids have not been systematically studied like dislocations in crystalline solids. Presently, we know little about their dynamics and their aggregated structures or emerging patterns/formations, much less about their relations to plasticity and ductile-to-brittle transition. For instance, recent researches have revealed that amorphous metals are capable of significant shear flow, which is mainly localized in extremely narrow shear bands that are about 10 nm in thickness. However, it has also been found that severe plastic instability may set in at the onset of plastic deformation, and corresponding localized heat in shear bands may lead to decohesion and subsequent crack propagation e.g. [17]. For instance, macroscale sized amorphous silica is often brittle at room temperature, despite of speculation of existence of minute plastic deformation or flow. 2

ACCEPTED MANUSCRIPT

EP

TE

D

M AN U

SC

RI PT

Molecular dynamics has been the main tool in modeling ductility or plasticity of amorphous materials e.g. plasticity in glass-forming materials with specific microstructure e.g. [18, 19, 20]. Most these studies are limited in a small region in nanometer scale. To capture plastic deformation in amorphous materials at macroscale, such as shear strain localization, without using continuum level empirical modeling is still a major technical challenge. In this work, we proposed a multiscale method that couples a coarse-grained Parrinello-Rahman (CG-PR) molecular mechanics with continuum finite element method, in which a representative sampling cell is used to derive macroscale stress, including flow stress, through the Cauchy-Born rule in each macroscale quadrature point in amorphous solids. The advantages of the proposed method are three-fold: (1) It is an atomistic-informed continuum approach that takes into account long range order flow stress; (2) the use of the representative sampling cell can take into account the effect of local disordered microstructure and precisely quantify its rearrangement and its defect evolution in a local cell of the amorphous solid, and (3) It is an atomistic interaction based method that does not have any empirical constitutive relation except that interatomic potential model may be semi-empirical.

AC C

Figure 1: Representative sampling cell with different molecular clusters (blue, green, yellow, orange, red, and silver) at different quadrature points in different elements.

The rests of the paper are organized into four sections: Section 2 provides detailed formulations of the proposed multiscale method namely, Multiscale Shear-TransformationZone Model. Section 3 focuses on discussion of numerical modeling, which presents several validation and prediction examples of plastic deformation in amorphous materials. In Section 4, we discuss some remaining issues. The continuum finite element formulations and the detailed molecular dynamics simulation process are provided in the Section 5. 2. Multiscale Shear-Transformation-Zone Model In this study we choose the Lennard-Jones binary glass (LJBG), as the amorphous material model to conduct the multiscale study, because it is a primary model for bulk metallic 3

ACCEPTED MANUSCRIPT

glasses [14]. The LJBG model is composed of two different sizes of atoms (A and B). LJBG model is a generic model of glassy materials, which has been extensively used to model glass formation of binary alloys, and it has been thoroughly studied to understand plastic deformation of metallic glasses e.g. [6, 21, 22].

TE

D

M AN U

SC

RI PT

2.1. Representative sampling cell (RS-cell) For the proposed multiscale amorphous model, a representative sampling cell (RS-cell) of 250 atoms of statistically representative configuration is embedded in every quadrature point in each element. In Fig. 1, we show a representative sampling cell for LJ-binary system, in which there are several atom clusters that are marked with different colors, i.e. blue, green, orange, yellow, pink, purple, red, and silver. An advantage of the multiscale method is its ability to capture microscale rearrangements of atom clusters, and subsequently translates these microscale molecular configuration change to the macroscale inelastic deformation.

EP

Figure 2: Illustration of kinematic mapping, configuration spaces of the coarse gained Parinnello-Rahman (CG-PR) method.

AC C

Therefore, a main function of the RS-cell is to host potential shear-transformation-zones in every macroscale finite element quadrature points. Before constructing such representative sampling cells, we have systematically conducted a series of molecular dynamics simulations (see Section 5.1), and we have found that the 250-atom RS-cell is the minimum size of the RS-cell that can reproduce mechanical responses of a larger atomic model over 30,000 atoms. In fact, a systematic study has been conducted to validate the RS-cell multiscale approach with the full-atom molecular dynamics approach for mechanical responses of macroscale amorphous solids under various applied loading conditions for several amorphous material models (see [23] for details). The main technique of the multiscale method is to create an RS-cell that can host potential shear-transformation-zones, and then embed the RS-cell into the representative macroscale material point, e.g. the site of a quadrature point in a macroscale finite element 4

RI PT

ACCEPTED MANUSCRIPT

SC

Figure 3: Schematic illustration of uniaxial tension and simple shear tests for a cubic specimen of binary amorphous solid.

AC C

EP

TE

D

M AN U

or finite element mesh. By doing so, during applied loadings, the amorphous material model can automatically produce plastic deformation due to the progression of the sheartransformation-zone and its manifestation at microscale. In specific, when an amorphous material undergoes deformation, each macroscale material point is stressed or each RS-cell is stressed, and when an RS-cell is stressed, the molecular clusters inside the RS-cell start to rearrange their relative positions, as shown in Figs. 1 and 6, in which a macroscale cubic specimen is under shear deformation. To quantify the shear-transformation-zone, we initially divide the RS-cell into several distinct regions, or clusters. To distinguish these clusters, we make the atoms in theses clusters in eight different colors, Blue (1), Green (2), Yellow (3), Orange (4), Pink (5), Purple (6), Red (7), and Silver (8) colors. There are several criteria to group different atoms into different groups. For example, dividing the RS-cell into 8 equal regions (volume), and all the atoms in a chosen region will be labeled in the same color. More relevantly, we used the k-means clustering algorithm [37] to arrange all the atoms into different clusters according to their affinity. The k-means clustering technique is an unsupervised machine learning algorithm that is extensively used in artificial intelligence and data-driven computations. This technique may potentially allow us to define and quantify defect initiation and evolution in amorphous solids, which will be reported in a subsequent paper. From Fig. 1, one can see that when the stress states in different elements are different, the relative positions, volumes, as well as the topological shapes of these clusters (colored atoms) are completely different. We also tend to believe that the void space among these different regions are also different or change under different stress states, which may be quantitatively computed by using the multiscale formulation described in the next Section. The shear-transformation-zone transition in amorphous solids is understood as a deformation process that is a collective motion of a group of atoms, which rearrange to achieve a local characteristic shear strain, γc , under an applied shear loading. A popular belief or hypothesis of the current STZ theory is that one may be able to capture STZ transition without the need of knowing the details of local rearrangement of individual atoms. However, the conventional STZ transition criteria is either based on the empirical measurement 5

ACCEPTED MANUSCRIPT

RI PT

or qualitative thermodynamics argument. The proposed multiscale shear-transformationzone approach is based on a multiscale molecular mechanics, which can automatically and precisely acquire a physical or “atomistic” STZ transition during computer simulations on the fly, according to the original concept in the physical-based STZ theory. The detailed technical procedure is explained below.

EP

TE

D

M AN U

SC

2.2. A coarse-grained Parrinello-Rahman (CG-PR) method In this section, we first introduce a coarse-grained Parrinello-Rahman (CG-PR) method, which was recently developed by the present authors for multiscale simulations of amorphous materials [23]. The CG-PR method has been carefully examined and validated by comparing its simulation results with those of lager sized models of MD simulations (see [23]). By using CG-PR method in conjunction with the RS-cell technique, we can simulate sheartransformation-zone transition, and the related plastic deformations in amorphous materials. In many multiscale finite element method for crystalline lattice materials, the CauchyBorn rule (CBR) is applied to an unit cell of a crystal, which is usually “inserted” into quadrature points within an element, and this procedure can be used to obtain macroscale stress at those quadrature points, e.g. [24, 25]. For amorphous materials, there is no periodically distributed unit cells. In the proposed CG-PR method, we employ a CBR procedure with applying a representative sampling cell (RS-cell) to extrapolate the stressstrain relation. Note that we cannot, however, directly apply CBR to derive its macroscale constitutive relation, since the amorphous material has no definite lattice structure. In order to construct a viable CBR procedure in amorphous materials, we have developed the following coarse-grained Parrinello-Rahman method, which originates from the well-known Parrinello-Rahman molecular dynamics [26, 27]. We first select the RS-cell of the amorphous material by identifying its basic microstructure and configuration as shown in Fig. 1. The RS-cell selected in multiscale simulation for the LJBG contains 250 atoms as mentioned above. We may define an RS-cell shape tensor as h(0) = [a, b, c], (1)

AC C

where a, b, c are three edge vectors of the unit cell at initial time, as suggested by Parrinello and Rahman [26]. Following [26], we may use the cell shape tensor h to re-scales the initial positions of atoms, i.e. Ri = h(0)Si , i = 1, 2, · · · Nc , (2) where Si is the scaled atom positions (see [27] for discussions). Considering the deformation mapping ϕ : R → r, one can define the deformation gradient as ∂ϕ F= . ∂R Applying the Cauchy-Born rule to approximate the deformation mapping ϕ, we can use the Cauchy-Born rule to write ri = FRi = F · h(0) · Si = h(t)Si , 6

(3)

ACCEPTED MANUSCRIPT

where h(t) = F · h(0) as shown in Fig. 2. In PR-MD [26], one solves the motion of Si (t) (or Ri (t)) for each atom i, and the evolution of cell shape h(t) (or F(t)) for each unit cell by using the following equations. Thus the equations of motion can rewritten as, X V 0 (rij )  j6=i

mi rij

˙ ·R ˙i (Ri − Rj ) − C−1 · C

RI PT

¨i = − R

¨ = −Ω0 F(S virial − S ext ) , WF

(4) (5)

SC

where W is the effective mass, or moment of inertia, of the cell, and S ext is the prescribed external second Piola-Kirchhoff stress; Ω0 is the volume of MD cell in the referential configuration, and  X V 0 (rij ) 1 X ˙i⊗R ˙i+ −mi R (Rj − Ri ) ⊗ (Rj − Ri ) Ω0 i rij j6=i  X V 0 (rij ) X 1 = −mi S˙ i ⊗ S˙ i + h0 · (Sj − Si ) ⊗ (Sj − Si ) · hT0 , Ω0 rij i j6=i

M AN U

S virial =

(6)

TE

D

where V (rij ) is the atomistic potential of the amorphous solid, and rij = |rj − ri |. For the detailed discussions and applications of the latest PR-MD, readers may consult [27, 28, 29]. For the coarse grained PR method (CG-PR), we do not solve each Si explicitly. On the other hand, we use multiscale finite element method to solve F(tn ) for each element at time tn = n∆t. In CG-PR computations, the atom configurations (rgi (tn )) in the current time step n is guessed from previous configuration by using the Cauchy-Born Rule (CBR), rgi (tn ) = F(tn ) · Ri (tn−1 ) = F(tn ) · h(0) · Si (tn−1 ),

(7)

AC C

EP

where, h(0) is the initial RS-cell shape tensor, Ri (tn−1 ) is atom coordinate in referential configuration at time step n − 1, and Si (tn−1 ) is normalized coordinates of atoms in the RS-cell. Once we have initial guess of coordinate as rge i (tn ), the configuration is optimized to minimize energy by using conjugate gradient method with periodic boundary condition in all direction. ropt i (tn ) = argminV (ri ) ,

(8)

with the starting guessed configuration rge i (tn ). In computations, we first use the CBR to predict the deformed configuration and stress, i.e. rge i = FRi , P(Xi ) =

and

ge 1 X 0 ge rij ⊗ Rij 1 X ∂V  = V (rij ) . ge 2Ω0 k ∂F 2Ω0 j rij

7

ACCEPTED MANUSCRIPT

To seek the optimal condition, ropt = argmin i

1 X V (rij ), 2Ω0 j

opt X ∂V (rij ) j

∂ropt i

=0 →

n X ∂V (rij ) j

∂rni

+

RI PT

we consider n X ∂ 2 V (rij ) j

∂rni ∂rni

· ∆rni ≈ 0 .

(9)

X ∂ 2 V (rn ) −1 X ∂V (rn )  ij ij =− , · n n n ∂ri ∂ri ∂ri j j

M AN U

∆rni

SC

By using Eq. (9), we can solve ∆rni , so that we can move to the next atom configuration by rn+1 = rni + ∆rni where i

opt with r0i = rge i . Once we have found ri , we can find

Sopt = (F · h)−1 · ropt i ,

and Ropt = hSopt . i i

(10)

Finally we can find the correct 1st Piola-Kirchhoff stress as follows, P=

1 ∂W 1 X 0 opt ropt ⊗ Ropt = V (rk ) k opt k , 2Ω0 ∂F 2Ω0 k rk

(11)

AC C

EP

TE

D

where k is the bond number index. In passing we note that the Cauchy-Born rule is essentially a procedure to obtain elastic stress state under an affine deformation. The CG-PR approach adopted here has an optimization part, which lets the system relax into an energy minimum configuration by optimally rearranging of molecular clusters’ position inside the RS-cell. Thus, it will produce inelastic stress responses at macroscale. It may be noted that we may be able to optimize cell shape tensor h by minimizing the potential energy V (r) in terms of h as well. In fact, the optimization procedure can be expressed as, 1 X V (rij ) . hopt = argmin 2Ω0 j If we use rn = Fhn Sni ,

opt X ∂V (rij ) j

∂hopt

=0 →

m X ∂V (rij ) j

∂hn

+

n X ∂ 2 V (rij ) j

∂hn ∂hn

· ∆hn ≈ 0 .

It is thus possible to obtain hn+1 = hn + ∆hn , and X ∂ 2 V (rn ) −1 X ∂V (rn )  ij ij ∆h = − · . n n n ∂h ∂h ∂h j j n

8

ACCEPTED MANUSCRIPT

RI PT

It is noted that this algorithm may not work, if we vary the cell shape tensor independent with the change of element shape, because the stress may diminish (relax) to reduce the force acting on atoms, and it then leads inconsistent density or volume change in each element. However, the change of det{h} may be taken as a “free volume” parameter as that of the internal variable based STZ theory of visco-plasticity to describe rheology of plastic flow e.g. [3, 18, 30], which will be discussed in a separated paper. In the present case, we restrict ourselves to optimize internal atom coordinate by using Eqs.(7)-(11) while maintaining the same cell shape tensor h = const., i.e. we do not consider the dilatation of the shear-transformation-zone [13].

SC

3. Multiscale Simulations

M AN U

Using the CG-PR method discussed in the above Section, we have conducted a series of simulations on plastic deformations of LJ binary glass material. In this Section, we discuss these simulation examples in details.

AC C

EP

TE

D

3.1. Example I: Uniaxial deformation In the first example, we use CG-PR method to simulate mechanical response or stressstrain response of LJBG under uniaxial deformation with both loading and unloading under different strain rates as shown in Fig. 3 (a). The dimension of the cubic specimen is: 100nm × 100nm × 100nm. There are total 100 tetrahedral elements used in multiscale finite element discretization. Each element contains only one quadrature point, and each quadrature point is associated with an RS-cell that contains 250 atoms. In the numerical simulations, we vary the loading/unloading strain rate in a range of ε˙ = ±2.0 × 107 s−1 ∼ 2.0 × 106 s−1 . Fig. 4 compares different loading and unloading stress-strain relations under different strain rates for the LJ binary amorphous solid, in which the cubic specimen is subjected to uniaxial deformation at different unloading strain levels, such as ±0.04, 0.06, 0.080.1, 0.15 and 0.2 strains, with the same strain rate. From Fig. 4, it is clear that lower strain rate decreases the yield stress in stretched deformation. Even at ε = 0.04 unloading strain, the stress-strain (S-S) relation curve is not the same as the original S-S curve for both tension and compressive deformations. This is because that there is not a definition yield stress after the first unloading for LJBG materials, and plastic deformation always occurs in the specimen as the S-S curves show work-softening less than ε = ±0.04. The discrepancy of between the S-S curve in reloading and the original S-S curve is larger when the specimen is re-deformed under larger strain. In addition, in Fig. 4, one can see that the residual strain i.e. the permanent strain measured when the stress becomes zero under high strain rate loading is larger than that under lower strain rate loading condition, indicating that slower deformation promotes plastic deformation more, because it has more time for atoms to rearrange. If one carefully exams Fig. 4, one may find that the stress-strain (S-S) relation in tension and compression is not completely symmetric, and there are noticeable differences in the yield stresses in different directions as well as the shape of the S-S curve in tension and 9

ACCEPTED MANUSCRIPT

D

M AN U

SC

RI PT

compression parts, which may be due to the Bauschinger effect in amorphous materials. According to the simulation results, it is fair to say that the CG-PR method enables us to simulate history-dependent plastic deformation of amorphous solids with different strain rates without assuming any empirical constitutive relations on continuum methods.

TE

Figure 4: Loading-unloading stress-strain loop under uniaxial tension at with different unloading strains and different strain rates obtained by using CG-PR method: (A) ε˙ = ±2.0 × 108 s−1 , and (B) ε˙ = ±2.0 × 107 s−1 (The dotted lines are unloading paths).

AC C

EP

In this work, in order to reflect and capture the randomness nature of amorphous materials, we have chosen different initial configurations for the RS-cells for different elements. To do so, three hundred different configurations were randomly selected (generated) by using MD simulations following the same melt-and-quench procedure to form a RS-cell database. In the multiscale simulation, one of these configurations is assigned to a quadrature point in an element of the cubic specimen, which is discretized by either 100 elements (100nm × 100nm × 100nm specimen), or 288 elements (100µm × 100µm × 100µm specimen), respectively. We compare the stress-strain relation of the specimen modeled with multiple different RS-cells with that modeled by using only one RS-cell but with randomly orientated initial configurations in the small specimen with the dimension 100nm × 100nm × 100nm. As shown in Fig. 5, the two stress-strain relations are almost the same. This indicates that it may be possible to use only one RS-cell with random different orientations to mimic amorphous random microstructure without losing physical fidelity, and this will save a lot of computation cost, because preparing different RS-cells for all different elements is time10

SC

RI PT

ACCEPTED MANUSCRIPT

M AN U

Figure 5: Uniaxial tention of (A) 100 nm size cubic discretized by 100 elements with a randomly rotated RScell, (B) 100 nm size cubic discretized by 100 elements with 100 RS-cells, (C) 100 µm size cubic discretized by 288 elements with 288 RS-cells. Loading strain rate for (A) and (B) is ε˙ = ±2.0 × 108 s−1 , while ε˙ = ±2.0 × 105 s−1 for (C).

AC C

EP

TE

D

consuming. Moreover, we also examined the effect of the specimen size and time scale by conducting the numerical simulation in a cubic specimen with side length L = 0.1mm by using 288 elements, which is 1000 times larger than the small cubic specimen. Again in this simulation, we randomly select an RS-cell from the RS-cell database of three-hundred different configurations, respectively. In this case, the loading strain rate is ε˙ = 2× 104 s−1 , which is three orders of magnitude slower than the previous two other cases. However, one may see from Fig. 5 that the overall profile of stress-strain relation under uniaxial stretch for the large specimen is similar or consistent with that of the previous two cases, which implies that we can simulate slower strain rate for larger size systems without any difficulties as long as deformation speed (L × ε) ˙ is within the same order of magnitude. More importantly, this result confirms that the proposed method is indeed a realistic multiscale method for continuum scale modeling. Note that here one may observe stress perturbation in the initial stage of the uniaxial deformation simulation. This may be attributed to the initial stress inside the RS-cell, because, in this method, we employ the conjugate gradient method to minimize energy without controlling stress of the RS-cell. Indeed, the initial stress fluctuation is not found in the case of shear deformation test, because the initial shear stress inside the RS-cell is negligible. Eliminating initial stress at energy minimum stage remains as one of the issues of this method. 3.2. Example II Shear deformation In the second example, we perform a simple shear test for the same cubic specimen under the following three loading and unloading strain rates: γ˙ = 1.0 × 108 s−1 , 1.0 × 107 s−1 , 1.0 × 106 s−1 . The same small cubic specimen used in the uniaxial tension test is used the simulation (see Fig. 3 (b)). In each element, a same RS-cell is embedded with a random (different) orientation, to mimic the random amorphous configuration for each quadrature point. 11

ACCEPTED MANUSCRIPT

TE

D

M AN U

SC

RI PT

We plot the inelastic stress-strain relations of LJBG material under pure shear loading/unloading with different strain rates in Figs. 6 (A), (B), and (C). From Fig. 6, one may find that under the shear loading condition, it is more obvious that the yield stress decreases when loading strain rate is small or slower, which is in a good agreement with that of MD simulations [23]. Again, the residual shear strain measured under the lower strain rate loading condition is larger than that under the higher strain rate loading condition. To confirm the validity of the multiscale shear-transformation-zone method, we display the relative positions, geometric shapes, and topological structures of all molecular clusters in a chosen RS-cell at different shear strains in Fig. 7. From Fig. 7 one can clearly see the spatial evolution of molecular clusters and the changes of their relative positions or arrangements corresponding to different shear strain levels or applied shear stress states. This clearly demonstrated that CG-RP method can automatically acquire shear-transformationzone transition condition, and directly incorporate it into the macroscale flow stress calculation while maintaining the elastic stress level without any empirical or phenomenological modelings, except for interatomic interaction model.

EP

Figure 6: Loading-unloading stress-strain loop under cyclic shear deformations at different unloading strains and under different strain rates by using CG-PR method: (A) γ˙ = 1.0 × 108 s−1 , (B) γ˙ = 1.0 × 107 s−1 , (C) γ˙ = 1.0 × 106 s−1 .

AC C

Another unique feature of the present method is that it can capture local deformation or local event at microscale owing to its possession of atomistic information inside each RS-cell. It is thus worth demonstrating that we can deduce local deformation state in the specimen by using the atom positions in RS-cells (without using finite elements). To quantify deformation state in the amorphous specimen, we first simulate shear deformation of the cubic specimen by employing multiple RS-cells, which are randomly selected from three hundred different configurations prepared in MD simulations, to all elements, individually. Then, local deformation gradient at each atom position Fatom in every RS-cells can be evaluated based on a. numerical procedure that is discussed in Section 5.3 and also in References [32, 27]. Fig. 8 (A) displays the distribution of F13 at atom position in the range of 0.05 to 0.3 macro shear strain when the cubic specimen is deformed with the stain rate γ˙ = 1 × 107 s−1 . When small strain is applied most of the atoms follow to macroscopic deformation, whereas the local deformation state varies widely under finite deformation conditions, indicating that essentially 12

SC

RI PT

ACCEPTED MANUSCRIPT

M AN U

Figure 7: The multiscale connection between microscale shear-transformation-zone evolution in an RS-cell and the macroscale stress-strain relations under pure shear loading. Change of relative positions, local topological arrangements, and configuration of different molecular clusters in an RS-cell corresponding to different strain levels: (a) the local region is centered in a same “blue atom”, and (b) the local region is centered in a same “yellow atom” in a real simulation.

EP

TE

D

amorphous materials do not obey affine deformation. Based on this result, we may expect that local deformation state in amorphous materials is fundamentally inhomogeneous and more localized under large macroscale deformation. In addition, interestingly, one may find significant differences in deformation distribution among different shear strain rates as seen in Fig. 8 (B). The distribution profile is almost symmetrical for the faster strain rate, but it is antisymmetric for the slowest case, γ˙ = 1 × 106 s−1 . This suggests that deformation is more localized in case of slower deformation, because it has more time to relax to an equilibrium state. These characterizations obtained by using the multiscale coupling method provide us insightful information on how to bridge the multiscale STZ modeling and the phenomenological STZ modeling, which will be examined further in subsequent studies.

AC C

3.3. Example III: Cyclic loading and cyclic plasticity A challenge problem in amorphous plasticity is to predict plastic damage evolution under cyclic loading. Despite of lack of emergence of defect patterns, plastic damage in amorphous can evolve under cyclic loading and provides a source of kinematically irreversible damage accumulation in amorphous materials, exhibiting fatigue failure under cyclic loading conditions well below the yield stress [19]. In this example, we used CG-PR method to calculate mechanical responses of LJBG specimen under cyclic loadings. Fig. 9 compares the stress-strain hysteresis loops under the uniaxial cyclic loadings with that under the shear cyclic loadings in two strain rate conditions: ε˙ = 2 × 107 s−1 and γ˙ = 1 × 106 s−1 . The stress-strain relations under cyclic shear loading conditions clearly show that the yield strength decreases. 13

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 8: Distribution of local deformation at each atom position Fatom during shear deformation. (A) 13 Distribution at different macroscale deformation with simple shear loading rate γ˙ = 1 × 107 s−1 . (B) Comparison of distributions at 0.2 shear strain among (a) γ˙ = 1 × 108 s−1 , (b) γ˙ = 1 × 107 s−1 and (c) γ˙ = 1 × 106 s−1 .

AC C

EP

TE

D

Moreover, in Fig. 9, one may find that the loading curves after the second cycle all become curly, which strongly implies the change of the yield stress as well. In fact, by examining Fig. 4, one may also find that both the yield stress values under tension and in compression loadings have changed depending on unloading strains. All these observations point to a possible Bauschinger effect for LJBG amorphous material, which is similar to the Bauschinger effect in amorphous materials observed before [36]. On the other hand, in metal plasticity, the Bauschinger effect refers to the anisotropic plasticity yielding. After carefully examining Fig. 4, one may find that both the tension and compression yield stress decrease simultaneously for LJBG, which is somewhat different from the Bauschinger effect in metallic materials, in which when tension yield stress increases while the yield stress in compression decrease. It seems to us that the yield surface for LJ binary glass is shrinking during each cyclic loading. That is: the yield stresses decrease for both tensile and compression loadings. To understand the atomistic origin of this phenomenon, we may need to exam the evolution and changes in all molecular clusters, which is beyond the scope of this work, and it will be reported in a separated paper. 3.4. Example IV: Shear band formations Shear banding is a major inelastic failure mechanism of bulk metallic glasses, e.g. [31]. To demonstrate the applicability of shear-transformation-zone theory, Rycroft and Gibou [33] have used the conventional shear-transformation-zone theory to calculate shear strain location formation, i.e. shear band formation, in an amorphous material [33]. To demonstrate the capacity of multiscale shear-transformation-zone model, we apply the multiscale 14

SC

RI PT

ACCEPTED MANUSCRIPT

(A)

(B)

M AN U

Figure 9: Cyclic loading/unloadling hysteresis loops in both uniaxial stretch and shear of the LJ binary glass. (A) Uniaxial stretch: loading rate ˙ = 2 × 107 s−1 , and (B) Simple shear: loading rate γ˙ = 1 × 107 s−1 .

TE

D

CG-PR method to simulate shear band formation in LJBG material. Moreover, in order to validate the multiscale shear-transformation-zone modeling, we conduct a comparison (controlled) molecular dynamics simulation of shear band formation in LJBG material at nanoscale, which is discussed in Section Methods. In the CG-PR simulation, following the example done by Rycroft and Gibou [33] , we conduct tension test on a macroscale shingle-notched specimen that has dimension 3mm × 0.1mm × 1mm. The initial notch shape has a width Lw = 0.8Lz and height Lh = 0.6Lz , and it is made created by removing the region, |X − X0 | 5

Lw cos(2Z/Lh − 1), 2π

(12)

AC C

EP

which is the same as the numerical model of Rycroft and Gibou [33]. Note that in Eq. (12) X0 is coordinate of center of the specimen in X direction, and Lz is width of the specimen in Z direction. To access the accuracy of the proposed method, we prepared two different FEM meshes, whose numbers of elements are 3,396 and 9,344, respectively, to investigate the convergence property of the multiscale method. In each element of a given FEM mesh, we randomly select one of 300 RS-cells, and embedded it into the element. Thus for the entire macroscale specimen, there are total 849,000 and 2,336,000 atoms used in the associated multiscale finite element method for the coarse and fine models, respectively. The discrete equation of motion of multiscale continuum is integrated with the time step 0.1 ns (see Section Methods). Uniaxial stretch boundary condition with 30 m/s of velocity, which corresponds to 2×105 s−1 strain rate, is applied to the upper and lower edges of the specimen. As shown in Figs. 10 (A), (B), the initial notch is gradually spreading, and the top side is bending. Consequently, V-shaped shear band is generated from the initial notch to the top side. At 0.08 strain, the shear band is clearly visible, and it is further enhanced at higher 15

ACCEPTED MANUSCRIPT

AC C

EP

TE

D

M AN U

SC

RI PT

strain 0.1. This result qualitatively agrees with nanoscale MD simulation results, which are displayed in Figs. 11 (A) and 11 (B). Comparing between MD and CG-PR results, it seems that the width of the shear band obtained by CG-PR multiscale method is broader than that obtained in MD simulations, because the coarse FEM mesh cannot resolve the resolution of the shear band. On the other hand it is almost identical if we use fine scale mesh model. Although the angle between the two shear band is slightly narrower in the case of CG-PR simulations than that in MD simulations, the results obtained by using CG-PR method clearly show localized effective shear strain in front of the notch. We can conclude that the CG-PR multiscale finite element method, which describes macroscale material response by using atomistic representative sampling cells, can capture shear band formation in glassy materials without using any empirical constitutive relation, during the plastic deformation process.

Figure 10: Shear band formation in the single-notched specimen of Lennard-Jones binary glass model by using CG-PR method. (A) Coarse scale model with 3,396 elements, (B) Fine scale mode with 9,344 elements,

4. Discussions In this work, we proposed a multiscale approach, i.e. CG-PR multiscale finite element method, to simulate plastic deformation of amorphous solids. By using CG-PR method, we have shown that during the mechanical loading the atom clusters inside each representative sampling cell can re-organize and rearrange their positions, which leads to permanent 16

ACCEPTED MANUSCRIPT

5. Methods

EP

TE

D

M AN U

SC

RI PT

plasticity deformation. To the authors’ best knowledge, this is the first successful multiscale simulation reported that uses the Cauchy-Born rule to capture inelastic deformation in amorphous solids, which is in contrast or in parallel with the similar simulations performed for the crystalline solids, e.g. [34]. In this study, we have observed the following facts: (1) The presented simulation results clearly demonstrated that the CG-PR multiscale finite element method can capture the plastic deformation of the binary amorphous solids at macroscale and room temperature. However, the simulation also shows that the simulated plasticity or ductility behavior depends on microscale atomistic potential, temperature, and specimen size. We have conducted similar numerical computations for amorphous Silica via CG-PR finite element method, and there is no significant plastic deformation in the simulation; (2) Each RS-cell may be viewed as an aggregated assemble of several loosely connected or cross-linked atom clusters, which we painted with different colors. We have shown that the topological arrangement of these atom clusters will change corresponding to different flow stress states, which indicates that the proposed CG-PR finite element may also serve as an intrinsic defect-mediate computational plasticity model, or defect dynamics model in amorphous solids, and (3) the change of the volume of RS-cell element i.e. det(δh) 6= 0 may be viewed as an internal variable in this multiscale shear-transformation-zone theory. As pointed out early, we may be able to develop a multiscale visco-plasticity theory based on the present multiscale CG-PR formulation by taking into account the effect of the δh tensor of the RS-cell. At the moment, this research project is still a work in progress. In this work, we first report the theoretical and computational formulations of the method and some of the early results. In fact, to completely capture the plastic deformation process inside the amorphous solids and the related ductile-to-brittle process, we need to understand temperature effect during the irreversible non-equilibrium process. A coupled thermal-mechanical CG-PR multiscale is also a work in progress, which will be reported in a subsequent paper.

AC C

5.1. Molecular dynamics simulation All the molecular dynamics simulations reported in this paper are conducted by using the open-source molecular dynamics code LAMMPS [35]. In the simulations reported in this work, the ratio of two atoms is 80:20 and the parameter sets for the Lennard-Jones potential are as follows; energy depths are εAA = 1.0ε0 , εBB = 1.5ε0 and εAB = 0.5ε0 . Those of particle size are σAA = 1.0r0 , σBB = 0.8r0 and σAB = 0.88r0 , respectively. The interatomic interaction is truncated at 2.245r0 . In order to connect the binary system with real materials, we presume that ε0 and σ0 are 2.4503 kJ/mol and 3.385 ˚ A, respectively. In this work, all RS-cells are prepared by using the following melt-and-quench process: the LJ binary glass system was initially configured in a dilute gas state to avoid particle overlap, and then the system was compressed by applying high pressure. Following the 17

SC

RI PT

ACCEPTED MANUSCRIPT

M AN U

Figure 11: Shear band formation in the single-notched specimen of Lennard-Jones binary glass model by using MD simulation with (A) T = 0.01, and (B) T = 0.2.

AC C

EP

TE

D

equilibration at normalized temperature (T = 1.0), the system was quenched until T = 0.02 with quenching rate 1.96 × 10−4 in normalized L-J reduced unit. After a NPT simulation at T = 0.02 with stress free condition, the system was eventually optimized by using the conjugate gradient (CG) method. Time integration interval for L-J binary glass was 1.0 × 10−3 in the LJ reduced unit. For mechanical responses at room temperature, RS-cells with three different numbers of atoms, 250, 1000, and 30,000 are tested in comparison with allatom molecular dynamics simulation [23], and we found that the RS-cell with 250 atoms with reasonable accuracy, while attaining computational efficiency. For the MD simulation of shear band formation in LJBG materials, the atomistic model for single-notched specimen was assembled as a 6 × 1 × 2 cells aggregation of an unit cell, which is composed of 30,000 atoms. The unit cell model is prepared by the melt-and-quench process as well, and the length of the side (Lx ) is 9.8 nm. Therefore the size of the singlenotched specimen is chosen as 59.1 × 9.8 × 19.7 nm3 . An initial notch is created the same way as that of CG-PR finite element model or that of macroscale model of Rycroft and Gibou [33]. The total number of atoms is the single-notched specimen is 352,817. The MD shear band simulations of the single-notched specimen were carried out with NVT ensemble by applying periodic boundary condition only for Y direction (thickness direction). Atoms locating within 1 nm from both edges are fixed and constant velocity in X direction is applied to them. The velocity is ±0.01 in L-J unit, and velocities of the atoms are set to zero in Y and Z directions. We study two temperature conditions, T = 0.01 and 0.2, to understand temperature effect on the fracture mechanism. All MD simulations were carried out by using LAMMPS [35] with 0.001 time step in L-J unit. In Fig. 9 (A) and (B), sequences of simulation snapshots for atomistic configurations at different temperatures are visually compared with that of CG-PR simulation. The deformation gradient in the specimen is evaluated as follows. First, we divide the specimen into 60 × 5 × 20 rectangular boxes, then atoms locating within a certain distance from each 18

ACCEPTED MANUSCRIPT

M AN U

SC

RI PT

edge is assigned to the node. During deformation simulation, the nodal displacement is assumed as average displacement of the atoms assigned to the node. Once we know the nodal displacement, it is possible to evaluate deformation gradient tensor F as the same manner with finite element method by B-matrix, which is derivative of shape functions. Then F13 is averaged over five layers along Y direction to depict it in the X-Z 2D plane. Detail of the analysis algorithm can be found in [20]. As a result, at both temperature conditions, LJBG specimen has clearly a pair of shear band formations, which are diagonal lines from the initial notch to the other side. Even at relatively small deformation conditions, shear band appears and the shear deformation is enhanced as the specimen is stretched further. The overall shear band structure and change of the notch shape both agree well with those calculated by the transformation zone theory. In the case of lower temperature, T = 0.01, the deformation of the specimen is antisymmetric, resulting generation of anisotropic shear band. It might be because initial antisymmetry generated in the specimen can not be relaxed at such low temperature. Indeed, at T = 0.2, a pair of the shear band is almost symmetric. However, except for the symmetric shape, shear band structure does not change much by temperature. 5.2. Finite element formulation At continuum scale, the amorphous solid is treated as an effective medium modeled by nonlinear continuum mechanics and finite element method. To formulate finite element formulation, we first define total Lagrangian of the continuum system as,   (13) L = K − Wint + Wext ,

EP

TE

D

where Wext , Wint and K are the external potential energy, the strain energy of the continuum and total kinetic energy, respectively. In Eq.(13), Z 1 ˙ ρu˙ · udV (14) K= 2 ΩZ W (F)dV, (15) Wint = Ω

where

AC C

where W is the strain energy density functional of strain, ρ and u˙ are the mass density and velocity field of the continuum, respectively. The Hamiltonian principle can then be derived in terms of displacement variation, Z t1   δK − (δW int + δW ext ) dt = 0, (16) t0

Z

t2

Z

t2

Z

˙ dt, δKdt = ρu˙ · δ udV t1 Ω Z h Z h i i ∂W δW int = : δF dV = P : δF dV, Ω ∂F Z Z Ω ¯ · δudS. δW ext = − b · δudV − T

(17)

t1

Ω

δΩt

19

(18) (19)

ACCEPTED MANUSCRIPT

¯ is the traction vector on the In Eq.(19), b is the body force inside the bulk medium, and T surface δΩt , respectively. Consequently, the Galerkin weak formulation can be formulated in terms of element summation as follows,  Z Z  Z    ne ne ne h h h h ¯ ¨ · δu + P : δF dV = A T · δu dS , (20) ρ0 u ρ0 b · δu dV + A A e=1

Ωe

e=1

Ωe

Γt

RI PT

e=1

uh (X) =

nX node i=1

SC

where, Ωe is the domain of the element e; Γt is the traction boundary of the system. The superscript h used in the above equation is the standard notation in computational mechanics to represent a discretized displacement field corresponding to FEM interpolation field. Using FEM interpolation functions, we can approximate the displacement field as follows, Ni (X)di ,

(21)

M AN U

where nnode is number of node in an element. According to Eqs.(20) and (21), the discrete equations of motion for FEM node displacements can be derived as, ¨ + f int (d) = f ext , Md

(22)

TE

D

where, M is the mass matrix, and f int , and f ext are force vectors from bulk elements, and external force, respectively. These quantities are defined as follows, Z ne M= A ρ0 NTe Ne dV (23) e=1 Ωe Z ne int f = A BTe Pe (d)dV (24) e=1 Ωe Z Z o ne n ext T ¯ e dS , NTe T f = A Ne Be dV + (25) e=1

∂Γt

Ωe

AC C

EP

where the superscript T denotes the transpose matrix; ρ0 is material density, and Be is the element strain-displacement matrix, " # ∂Ne Be := . (26) ∂X In multiscale CG-PR method, or the multiscale shear-transformation-zone method, we integrate the macroscale dynamic equation (22) to find material deformation states. 5.3. Local deformation state analysis Local deformation state at each atom site was evaluated as the same approach in the particle mesh method [32] or multiscale MD method [27]. First, the moment matrix M of each atom i in initial configuration is calculated by using the initial coordination R as, Mi =

N X

ω(|Rij |)Rij ⊗ Rij

j=1

20

(27)

ACCEPTED MANUSCRIPT

RI PT

where, Rij = Rj − Ri and ω is a localized positive window function. A possible choice for the window function is the Gaussian function, ! 1 −x · x (28) ω(x) = exp (πh2 )3/2 h2

Ni =

N X

SC

where, h is the radius of the support and we set it to 5.0 ˚ A. To construct the microscale deformation gradient tensor at each atom position at time t (ri (t)), we need to evaluate a two-point second order tensor constructed by the following discrete sum, ω(|Rij |)rij (t) ⊗ Rij .

j=1

(29)

M AN U

According to Eqs. (27) and (29), a nanoscale deformation gradient Fi may be eventually estimated as Fi = Ni · M−1 i .

Acknowledgment

TE

D

The authors would like to thank Professor Takahiro Murashima of Tohoku University for his help using and integrating LAPPMS into the multiscale code. The authors are very grateful to Mrs. Mari Namba for her help to make the FEM mesh model.

AC C

EP

[1] F. Spaepen, A microscopic mechanism for steady state inhomogeneous flow in metallic glasses, Acta metallurgica, (1977), 25(4), pp.407-415. [2] A. S. Argon, Plastic deformation in metallic glasses, Acta metallurgica, (1979) 27(1), pp.47-58. [3] A.S. Argon and L.T. Shi, Development of visco-plastic deformation in metallic glasses, Acta Metallurgica, (1983), 31(4), pp.499-507. [4] E. R. Homer and C. A. Schuh, Mesoscale modeling of amorphous metals by shear transformation zone dynamics, Acta Materialia, (2009), 57, 2823-2833. [5] M. L. Falk and J. S. Langer, Dynamics of viscoplastic deformation in amorphous solids, Phys. Rev. E, 57, (1998) pp. 7192. [6] M.L. Falk, Molecular-dynamics study of ductile and brittle fracture in model noncrystalline solids, Physical Review B, (1999), 60(10), pp. 7062. [7] M.L. Falk and J. S. Langer, Deformation and failure of amorphous, solidlike materials, Annu. Rev. Condens. Matter Phys., (2011), 2(1), pp. 353-373. [8] J. S. Langer and L. Pechenik, Dynamics of shear-transformation zones in amorphous plasticity: Energetic constraints in a minimal theory, Phys. Rev. E, (2003), 68(6), No. 061507. [9] J. S. Langer, Dynamics of shear-transformation zones in amorphous plasticity: Formulation in terms of an effective disorder temperature, Phys. Rev. E, (2004), 70, No. 041502. [10] J. S. Langer, Shear-transformation-zone theory of plastic deformation near the glass transition, Physical Review E, (2008) 77(2), pp.021502. [11] C. M. Lee, K.W. Park, B.J. Lee, Y. Shibutani, and J.C. Lee, Structural disordering of amorphous alloys: A molecular dynamics analysis, Scripta Materialia, (2009), 61(9), 911-914.

21

ACCEPTED MANUSCRIPT

AC C

EP

TE

D

M AN U

SC

RI PT

[12] J. Luo, J. Wang, , E. Bitzek, J.Y. Huang, H. Zheng, L. Tong, Q. Yang, J. Li, and S. X. Mao, Sizedependent brittle-to-ductile transition in silica glass nanofibers, Nano letters, (2015) 16(1), pp.105-113. [13] M.Q. Jiang, G. Wilde, F. Jiang, and L. H. Dai, Understanding ductile-to-brittle transition of metallic glasses from shear transformation zone dilatation, Theoretical and Applied Mechanics Letters, (2015) 5(5), pp.200-204. [14] K. Zhang, H. Wang, S. Papanikolaou, Y. Liu, J. Schroers, M.D. Shattuck, and C. S. O’Hern, Computational studies of the glass-forming ability of model bulk metallic glasses, The Journal of Chemical Physics, (2013) 139(12), p.124503. [15] C. Zhong, H. Zhang, Q.P. Cao, X. D. Wang, D. X. Zhang, U. Ramamurty, and J. Z. Jiang, J.Z., Deformation behavior of metallic glasses with shear band like atomic structure: a molecular dynamics study, Scientific Reports, (2016) 6, p.30935. [16] H. B. Ke, J.F. Zeng, C.T. Liu, and Y. Yang, Y., Structure heterogeneity in metallic glass: Modeling and experiment, Journal of Materials Science and Technology, (2014) 30(6), pp. 560-565. [17] J.J. Lewandowski and A.L. Greer, Temperature rise at shear bands in metallic glasses, Nature Mat., (2008), 5, pp. 15-18. [18] Y. Qi, T. Cagin, Y. Kimura, and W. A. Goddard III, Molecular-dynamics simulations of glass formation and crystallization in binary liquid metals: Cu-Ag and Cu-Ni, Physical Review B, (1999), 59, pp. 35373533. [19] D. Chuang and C. A. Schuh, Atomistic mechanisms of cyclic hardening in metallic glass , Applied Physics Letters, (2012), 100, No. 251909 [20] S. Urata and Y. Sato, A study on the plasticity of soda-lime silica glass via molecular dynamics simulations, The Journal of Chemical Physics, (2017), 147(17), No. 174501. [21] C.D. Lorenz and M.J. Stevens, Fracture behavior of Lennard-Jones glasses, Physical Review E, (2003), 68(2), No. 021802. [22] F. Varnik, L. Bocquet, and J. L. Barrat, J.L., A study of the static yield stress in a binary Lennard-Jones glass, The Journal of chemical physics, (2004) 120(6), pp.2788-2801. [23] S. Urata and S. Li, A multiscale model for amorphous materials, Computational Materials Science, (2017), 135, pp. 64–77. [24] S. Li, and G. Wang, Introduction to micromechanics and nanomechanics, (2008) World Scientific Publishing Co Inc. [25] S. Urata and S. Li, Higher order Cauchy-Born rule based multiscale cohesive zone model and prediction of fracture toughness of Silicon thin films, International Journal of Fracture, (2017), 203(1), pp.159-181. [26] M. Parrinello, and A. Rahman, Polymorphic transitions in single crystals: A new molecular dynamics method, Journal of Applied physics, (1981), 52(12), pp. 7182-7190. [27] S. Li and S. Urata, An atomistic-to-continuum molecular dynamics: Theory, algorithm, and applications, Computer Methods in Applied Mechanics and Engineering, 306, (2016) 306, pp. 452-475. [28] P. Podio-Guidugli, On (Andersen)-Parrinello-Rahman molecular dynamics, the related metadynamics, and the Use of the Cauchy-Born rule, Journal of Elasticity, (2010), 100, 145-153. [29] S. Li and Q. Tong, A concurrent multiscale micromorphic molecular dynamics, Journal of Applied Physics, (2015), 117, No. 154303. [30] Y. Shi, M.B. Katz, H. Li, and M.L. Falk, Evaluation of the disorder temperature and free-volume formalisms via simulations of shear banding in amorphous solids, Physical review letters, (2007) 98(18), pp.185505. [31] F. Shimizu, S. Ogata, and J. Li, Theory of shear banding in metallic glasses and molecular dynamics calculations, Materuials Transactions, (2007), 48(11), pp. 2923-2927. [32] S. Li and W. K. Liu, Reproducing kernel hierarchical partition of unity, part I: Formulation and theory, International Journal for Numerical Methods in Engineering, (1999), 45(3), pp. 251-288. [33] C. H. Rycroft and F. Gibou, Simulations of a stretching bar using a plasticity model from the shear transformation zone theory, Journal of Computational Physics, (2012), 231, pp. 2155-2179. [34] D. Lyu and S. Li, Multiscale crystal defect dynamics: A coarse-grained lattice defect model based on crystal microstructure, Journal of the Mechanics and Physics of Solids, (2017), 107, pp. 379-410.

22

ACCEPTED MANUSCRIPT

AC C

EP

TE

D

M AN U

SC

RI PT

[35] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, Journal of Computational Physics, (1995), 117, pp. 1-19. [36] S. Karmakar, E. Lerner, and I. Procaccia, 2010. Plasticity-induced anisotropy in amorphous solids: The Bauschinger effect, Physical Review E, (2010) 82(2), p.026104. [37] T. Kanungo, D.M. Mount, M.S. Netanyahu, C.D. Piatko, R. Silverman, and A. Y. Wu, An efficient k-means clustering algorithm: Analysis and implementation, IEEE transactions on pattern analysis and machine intelligence, (2002) 24(7), pp. 881-892.

23

ACCEPTED MANUSCRIPT

A Multiscale Shear-Transformation -Zone Model and Simulation of Plasticity in Amorphous Materials

RI PT

Graphic Abstract

AC C

EP

TE D

M AN U

SC

1. We embed possible shear-transformation zone into the representative sampling cell (RS-cell) of LJ binary materials:

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

2. By using the proposed multiscale STZ model, we can predict inelastic deformation or stressstrain relation at continuum scale

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

3. By using the multiscale STZ method, we can predict cyclic plasticity,

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

4. We have successfully simulated shear band formation (strain localization) at macroscale (mm scale) by using multiscale STZ method: