Modeling and simulation of size effects in metallic glasses with a nonlocal continuum mechanics theory

1 Introduction

Metallic glass is a solid, amorphous material that was first produced in 1960 [1]. Most metallic glasses are alloys and often manufactured with the help of very rapid cooling to avoid crystallization. At first, scientists suggested that metallic glasses have a chaotic structure. However, lately, it is believed that they have some kind of semi-ordered structure [2].

Several types of metallic glasses with different distinctive properties exist. Zhang et al. [3] report on a metallic glass with an extremely low glass transition temperature and a polymer-like thermoplastic behavior. On the one hand, metallic glasses are stronger than their crystalline counterparts due to the lack of dislocations. On the other hand, this causes the metallic glasses to be more brittle.

However, there also exist exceptional types with a large ductility [4]. Further, metallic glasses are tougher than ceramics and have a greater elasticity, wear resistance, and corrosion resistance [2, 5]. In addition, their thermal and electrical conductivity is lower compared with crystalline materials.

Nowadays, metallic glasses are used or aimed to be used in sports equipment (e.g., golf club heads, golf balls, skis, baseball hats, tennis rackets), watches, medical devises (screws, pins, or plates for implantation into bones, scalpels), and automobile industry (because it is twice as strong as steel but lighter) (see also [6] for a list of possible applications). At the moment, production costs are still rather high, which limits the actual use. The mechanical properties of metallic glasses have been analyzed intensively, and consequently, there are several experimental investigations of metallic glasses [4, 7–15]. An extensive overview over the mechanical properties of bulk metallic glasses is given in [16].

Like crystalline metals, metallic glasses exhibit localized, i.e., heterogeneous, deformations. During deformation, shear bands form inside the metallic glasses. Following the explanation by Schuh and Lund [17], the shear transformation zone is the fundamental unit of plasticity in metallic glasses. In the shear transformation zone, small clusters of randomly close-packed atoms spontaneously and cooperatively reorganize under the applied shear strain. Self-assembly is the process that leads to the continued propagation of shear strains. This effect produces localized distortion within the material, which results in free-volume generation and thermal softening [16]. This triggers the formation of large planar bands of shear transformation zones, representing distinct shear bands in the material.

In the last years, there have been numerous experimental studies regarding the size effect in small-sized metallic glass samples with respect to the change in ductility, strength, or hardness (see, e.g., [18–24]). One major outcome of these studies is that the ductility increases with decreasing sample size. The main reason for this behavior is a delay or even suppression of shear localization with decreasing sample size [25]. Guo et al. [26] observed a stable growth of shear localization in their experiments with small sample sizes. This leads to much more plastic deformation before fracture, as observed in bulk metallic glasses, where usually catastrophic shear localization is obtained, see, e.g., [27]. Similar observations have been done in nanostructured materials [28–30]. Due to the small size of the specimen, the experimental preparation is very difficult and may affect the experimental results, such as processing-induced size effects [31]. Next to intrinsic effects, extrinsic effects (e.g., [21]) may also play a significant role. Therefore, modeling can help to further understand the mechanism and isolate the particular effects.

Due to the dominant failure mechanism of shear localization, which leads to a strain-softening behavior in metallic glasses, one encounters several problems in modeling. The application of classical continuum models result in loss of ellipticity of the boundary value problem, and a strong mesh dependency is observed. Several authors have dealt with these problems and obtained mesh-independent formulations (e.g., [32–34]), especially by the addition of gradient terms that preserve the ellipticity in the softening regime. In the works by Aifantis [35, 36], it was shown that the Laplacian in the plastic strain or a diffusive-like term in the internal variables settle the issue of shear band thickness and mesh-size dependence in the modeling of shear band formation.

For years, modeling the size-dependence of (poly-)crystalline materials, especially metals, has been of high research interest (see, e.g., [37–44]). Metallic glasses, which are amorphous metals and therefore noncrystalline, however, have hardly been studied in the framework of continuum mechanical approaches. Argon [45] present one of the first models of the plastic deformation in metallic glasses based on the introduction of the notion of the free volume. Anand and Su [46] present a finite-deformation, elastic-viscoplastic theory for metallic glasses. Purely theoretical approaches are presented by Vaks [47] and Huang et al. [48]. Steif et al. [49], Gao [50], and Thamburaja [51] computationally study the length-scale effects on the shear localization process in metallic glasses. A Ginzburg-Landau type of theory is introduced by Zheng and Li [52]. A rate-dependent theory for bulk metallic glasses, which is experimentally validated by the metallic glass named Vitreloy 1, is studied by Yang et al. [53]. The approach of [14] is based on the idea that metallic glasses have a similar behavior to granular materials. They assume that the motion of the atoms is comparable to the sliding of granules.

The purpose of this contribution is the modeling and analysis of size effects in metallic glasses. For this purpose, a nonlocal material model is proposed. To the authors’ knowledge, this is the first fully nonlocal continuum-mechanics-based approach attempting to model size effects in metallic glasses. Its mathematical formulation is presented in Section 2, and its finite element discretization in Section 3. A numerical investigation of the material behavior with respect to size effect, rate dependence, tension-compression asymmetry, and the influence of initial conditions is given in Section 4. The work ends with a brief summary in Section 5.

2 Mathematical model

The purpose of this section is to present the mathematical formulation for the modeling of metallic glasses. This is formulated in the framework of continuum thermodynamics (e.g., [54]) and rate variational methods (e.g., [44, 55]). To this end, let

The metallic glass consists of atoms of different sizes. This leads to a free volume¹ξ inside the material, which determines the inelastic deformation. The free-volume generation ξ is either induced by plastic shearing (i.e., plastic strain γ) and other mechanisms (i.e., diffusion, hydrostatic pressure, or structural relaxation), which are accounted for by ξ_m. For simplicity, the current formulation is restricted to quasi-static, infinitesimal deformation, supply-free, and isothermal processes.

Energy storage is represented by the general form ψ=ψ(E, E_p, γ, ξ, ∇ξ) of the free energy density depending in general on the strain tensor for small strains²E=sym ∇u the inelastic strain tensor E_p, and the internal variables γ, ξ, and ∇ξ.

In this contribution, the free energy density ψ is additively decomposed into three parts, i.e.,

The elastic part ψ_e is determined by

where

The plastic strain tensor E_p is decomposed into a spherical and deviatoric part

Following Thamburaja and Ekambaram [56] and Thamburaja [51], the deviatoric part is determined by

where N denotes the traceless unit vector of the plastic flow direction. This part describes the plastic deformation due to shear-like motion of localized atom groups [51].

The spherical part, which describes the plastic deformation due to free-volume generation, is determined by the total free-volume generation

and is not affected by an initial reference free volume ξ₀. Consequently, the plastic strain tensor is given by

Based on the modeling approaches of Demetriou and Johnson [57] and Heggen et al. [58], the total free-volume generation is decomposed as

where

The defect free energy density is determined by

(see also Thamburaja [51]), where a_def is a material constant representing the defect-free energy coefficient and ξ₀ denotes the fully annealed free volume.

The gradient part of the free energy density is assumed to be a quadratic function of the free-volume gradient ∇ξ:

depending on the fracture surface energy a_sf and the material length l.

Next, we are turning to dissipative kinetic processes, which are represented by the dissipation potential χ. The dissipation potential χ is additively decomposed into a part related to plastic strain γ and a part for all further mechanisms for free-volume generation

Power law relations are assumed for both processes,

where

Following the work of Thamburaja [51], the cohesion is determined by the evolution equation:

where c(0)=c₀ is the initial value of resistance, b is a unitless fitting constant, which is relevant for the rate of strain softening due to free-volume generation, and f₁=f₀ is a characteristic frequency. To model softening, the fitting constant b has to fulfill

The dissipation principle [54] is satisfied sufficiently in the bulk when the dissipation potential χ is nonnegative and convex in its nonequilibrium arguments (see, e.g., [55]).

With these basic constitutive relations, we now carry out the continuum thermodynamic variation formulation of the evolution-field relations of the model following [55]. The formulation begins with the following rate functional

The surface rate potential r_s:=ζ_s+χ_s consists of an energetic ζ_s and a dissipative kinetic χ_s part that are linear and nonlinear, respectively, in terms of the rates

The volumetric rate potential r_v:=ζ_v+χ_v is determined by the dissipation potential χ and the volumetric energy storage rate density ζ_v=ζ. The latter is determined by the free energy density ψ, which is written as

The energy storage rate density ζ is linear in the rates

The rate functional R is stationary with respect to all admissible variations of

which holds in the bulk and at the flux part of the boundary

The stationary condition of the rate functional R with respect to admissible variations of

which can be interpreted as generalized flow rules in the respective regions.

The stationary condition of the rate functional R with respect to admissible variations of

which can be physically interpreted as generalized free-volume generation rules in the respective regions.

Thus, the following formulations represent the strong form of the field equations in the bulk yielding

for the quasi-static momentum balance, the flow rule, and the free-volume generation rule, respectively.

Turning back to the energy storage rate density ζ, the form is written as

Here, we use the following definitions

where dev(σ) denotes the deviator of the Cauchy stress tensor σ, which is calculated as usual

Substituting Eqs. (11) and (20) into Eq. (19) results in the strong forms

The term

The evolution equation for the total free-volume generation ξ is obtained by inserting Eq. (8) into the third equation of Eq. (23):

3 Finite element discretization

Based on these relations, we formulate the finite-element algorithm for the initial boundary value problem to solve for the unknowns u, γ, and ξ. Thamburaja [51] uses a theory which provides the information on the integration points by a special numbering of the integration points where he is using the finite-difference method to solve for the Laplacian of the free volume ξ in the evolution equation of the plastic strain γ [second equation of Eq.(23)] and the free volume generation [Eq. (24)]. However, the use of this theory leads to problems when changing the numbering or the element sizes within the structure. For this reason a fully non-local method is proposed in this work.

We apply a dual mixed finite element solution method [62]: we locally solve for the plastic strain γ and the free volume ξ. At the global level, we solve for the displacement u and the free-volume generation gradient

Consequently, we work with the weak forms of the first equation of Eq. (23) and Eq. (25), i.e.,

satisfying the boundary conditions for the variations

as well as the second equation of Eq. (23) for the plastic strain γ, which we solve in its dual form

Domain

4 Numerical experiments

As opposed to crystalline metals, metallic glasses do not melt at a certain temperature but change reversibly from a solid to a liquid state at the glass transition temperature. This distinctive property is of great advantage because it helps to manufacture complex shapes out of the metallic glass. The material parameters used in this work are taken from the literature to reproduce the glassy state response of metallic glasses at low homologous temperatures. Homologous temperature describes the temperature as a fraction of its melting point temperature. As reported in [20], the plastic zone size (also referred to as fracture process zone f) lies between 100 nm and 100 μm. Samples smaller than the plastic zone will not exhibit brittle failure by catastrophic shear localization because of the stable propagation of shear bands. Therefore, shear bands become stabilized when the sample size is comparable to the plastic zone. However, the absolute value of sample size and material length scale is not of importance – rather its

The defect-free energy coefficient a_def as well as the resistance to free-volume generation due to mechanisms other than plastic strain (shearing) a_res are taken from [64]. The internal strain rate

Table 1

Material parameter values adopted for the mathematical model.

Parameter	Symbol	Value	Unit
Young’s modulus	E	100	GPa
Poisson’s ratio	v	0.4	–
Free-volume creation factor due to plastic strain	sγ	0.02	–
Fracture surface energy	asf	1	J m-2
Material length scale	l	10	μm
Defect-free energy coefficient	adef	3500	GJ m-3
Resistance to free-volume generation	ares	320	GJ m-3
Fully annealed free (reference) volume (at θ0=300K)	ξ0	0.06	%
Initial cohesion	c0	750	MPa
Material strain rate		0.001	s-1
Fitting constant	b	-300	–
Fitting constant	ϕ	0.2	–
Strain rate sensitivity	n	50	–
External applied strain rate		0.001	s-1

The value of fracture surface energy a_sf is taken from [65], the ones of the defect-free energy coefficient a_def and resistance to free-volume generation due to mechanisms other than plastic strain a_res from [64], and fully annealed free (reference) volume (at θ₀=300 K) ξ₀ from [53].

The numerical investigation is performed for two samples. The setup of the first sample S₁ is shown in Figure 1. Here, a rectangular with a size ratio of 2:1 is chosen, which is discretized by 40×80 bilinear plane-strain elements. The bottom face is fixed in the y (22)-direction as well as the node in the bottom left corner in x (11)- and y (22)-directions. The deformation is applied at the top face, which leads either to compression or tension depending on the applied loading direction. For the second field, micro-free boundary conditions h·n=0 (see [66]) are applied. A slightly lower value of the initial cohesion c₀ in the model represents a defect in the material and acts as nucleation for shear localization. If not stated otherwise, one element with a 1% lower initial cohesion c₀ is placed near the center of the samples.

Figure 1

Sample S₁ with dimension l₀×2l₀. The mesh consists of 40×80 elements. A defect is placed near the center to trigger the shear band in the simulation.

Sample S₂ is investigated as a second example, whose setup is shown in Figure 2. S₂ is a square discretized by 40×40 bilinear plane-strain elements and subject to plane-strain tension loading assuming periodic deformation and periodic boundary conditions for the free volume ξ. To this end, the bottom left corner of the square is fixed in the x (11)- and y (22)-directions, the top left corner in x-direction, and the bottom right corner in y-direction. A displacement rate

Figure 2

Sample S₂ with dimension l₀×l₀. The mesh consists of 40×40 elements. If not stated otherwise, a defect is placed near the center to trigger the shear band in the simulation. For this sample, periodic boundary conditions are assumed for the deformation and free volume ξ.

4.1 General model behavior

In Figure 3, the average stress-strain behavior of the model under compressive loading is shown for sample S₁. The applied boundary conditions lead to the fact that the average of all stress components (except σ₂₂) are zero, i.e., ∫_Bσ₁₁ dV≈0 and ∫_Bσ₁₂ dV≈0. The general model behavior is summarized as follows: First, the specimen behaves elastically. Next, the model predicts a macroscopically homogeneous plastic deformation that is characterized by a stress plateau (with slightly negative slope).

Figure 3

Sample S₁, general model behavior: average stress component |σ₂₂| over applied strain |E₂₂| under compressive loading.

This is followed by a stress drop in the simulation, which marks the start of shear localization, resulting in the formation of distinct shear bands. A plateau is reached once one shear band is fully developed. This general behavior is similar to results obtained in experiments [70] and molecular dynamic simulations [71, 72] for metallic glasses loaded under compression.

The spatial distribution of the free volume ξ at three characteristic applied compression strains (2%, 5%, and 8%) is displayed in Figure 4. In the elastic region (|E₂₂|<1.5%), no change in free volume ξ, plastic slip γ, or cohesion c is observed. After yielding, the free volume increases and shows a cross-structure, triggered by the lower initial cohesion c₀ near the center of the sample. After the stress drop (|E₂₂|≈3%), shear localization starts, the formation of one distinct shear band is noticed. The cross-structure, which was observed previously, vanishes on the cost of the formation of one distinct shear band. At |E₂₂|≈6%, the shear band is fully formed. With increasing deformation, the width of the shear band increases and the upper sample part slides along the shear band. However, no increase in the maximum value of free volume ξ is observed anymore.

Figure 4

Sample S₁, general model behavior: spatial distribution of free volume ξ at three applied compression strains |E₂₂|.

Figure 5 shows the spatial distribution of the cohesion c at the three characteristic applied compression strains (2%, 5%, and 8%). With increasing deformation and free volume ξ, the cohesion c is decreasing. The cohesion c describes the internal resistance of the atomic structure to sliding, therefore, a decreasing cohesion c promotes the development of shear bands. Cross-structured shear bands develop, which end up in one distinct shear band. As for ξ, the width of the shear band increases with increasing deformation and the minimum cohesion c is nearly constant between |E₂₂|=5% and |E₂₂|=8%.

Figure 5

Sample S₁, general model behavior: spatial distribution of cohesion c at three applied compression strains |E₂₂|.

The plastic strain distribution is shown in Figure 6 at three strain states. The plastic strain γ and free volume ξ are strongly coupled [see Eq. (24)]. Consequently, their behavior are similar. However, for |E₂₂|=8%, differences that are due to the gradient h are observed. The plastic strain distribution shows a stronger bulge distribution compared with the free-volume distribution.

Figure 6

Sample S₁, general model behavior: spatial distribution of plastic strain γ at three applied compression strains |E₂₂|.

To confirm the previous assumptions, in Figure 7, the development of the free volume ξ for different local elements of the sample are plotted. Four characteristic elements, two within the shear band (P₁, P₂), one within the side localized region (P₃), and one outside of any localized region (P₄), are chosen. The locations are shown in the spatial free-volume ξ distribution at |E₂₂|=8% of the meshed sample S₁. The side localized region refers to the region of the cross-structure at the beginning of the deformation (see Figure 4), which is not part of any distinct shear band.

Figure 7

Sample S₁, general model behavior: free volume ξ over applied strain |E₂₂| under compression loading for different elements of the mesh that shows the free-volume distribution ξ at |E₂₂|=8%.

P₁, element within shear band and with initial defect; P₂, element within shear band; P₃, element within side localized region; P₄, element outside localized regions.

After yielding, the free volume ξ increases for all elements in the same way. At |E₂₂|≈2.5%, differences in the free volume ξ between the elements are observed. The free volume in element P₄ already reaches a saturation value at this strain, whereas for element P₃, the value reaches a plateau at |E₂₂|≈3%. This strain state corresponds to the stress drop and the beginning of shear localization. For the development of the free volume ξ, nearly no difference is seen for the elements P₁ and P₂. The free volume reaches a constant value of ξ≈0.7% at |E₂₂|≈6%. Note that this is an outcome of the simulation and not due to any imposed saturation value. The observation of a saturation value for the free volume ξ agrees with assumptions and observations in the literature [46, 51].

Next, the development of the cohesion c is investigated (see Figure 8). The cohesion c shows the opposite behavior in comparison to the free volume ξ. With increasing applied strain |E₂₂|, the cohesion decreases. Analogous to the free volume ξ, the saturation values of c are reached after a certain amount of applied strains. Here, it is interesting to note the differences between elements P₁ and P₂. The initial cohesion c₀ of element P₁ is slightly smaller than in the rest of the structure. The development for elements P₁ and P₂ are similar up to |E₂₂|≈2.5%. Subsequently the cohesion c decreases for element P₂ until it nearly reaches the one from P₁. A saturation value of 100 MPa is reached for both elements.

Figure 8

Sample S₁, general model behavior: cohesion c over applied strain |E₂₂| under compression loading for different elements of the mesh that shows the cohesion distribution c at |E₂₂|=8%.

P₁, element within shear band and with initial defect; P₂, element within shear band; P₃, element within side localized region; P₄, element outside localized regions.

The presented nonlocal model for metallic glasses is able to model the stable growth of shear localization in submicron samples, resulting in a distinct shear band formation. The shear localization process starts at the specimen’s initial defect. With ongoing deformation, the saturation values of the free volume and cohesion are reached within the distinct shear band.

4.2 Length scale-dependent behavior

To investigate the length scale-dependent behavior, the

ratio of internal length l and external specimen size³l₀ is varied between 100 and 5000. The average stress-strain behavior is shown in Figure 9. Differences are observed in the strain range |E₂₂|≈3–6%, the so-called knee region [51]. An increased

ratio of internal length l to external specimen size l₀ can be interpreted as decreased sample size or increased fracture process zone size f. With increasing

ratio, the knee region diffuses, which implies that the mechanism of shear localization changes. Consequently, the shear localization process delays with an increase in

ratio.

Figure 9

Sample S₁, length scale-dependent behavior: average stress component |σ₂₂| over applied strain |E₂₂| under compression loading for different

ratios of internal length l and external specimen size l₀.

Figure 10 shows the spatial distribution for the different

ratios of the free volume ξ at 8% strain. For the smallest

ratio (=100), the most pronounced shear band is obtained, which represents a stable shear localization process. An increase in

leads to less pronounced shear bands and for

even no shear band is obtained. However, in this case, the formation of a bulge around the middle axis is observed. A further increase leads to a homogeneous deformation process (shear localization is suppressed) that is correlated to the size of the shear band nucleus. As reported in [73], the size of the shear band nucleus d for metallic glasses is around 20 nm. Therefore, a mechanism change occurs once

reaches the critical ratio

of the fracture process zone size f and shear band nucleus d, implying that the sample size is of the same size as the shear band nucleus and no shear band is nucleated in the sample. As shown in [51], for

unstable (catastrophic) shear localization is obtained. A stable shear band forms for

The main mechanism that is responsible for the observed size effect is the nonlocal back stress [cf. Eq. (24)], which is energetic. As noted in [51], the back stress is the main mechanism. It is physically sensible because the free-volume diffusivity is very low in an annealed, and therefore relaxed, metallic glass sample (for temperatures below glass transition temperature). The back stress leads to a strengthening of the shear band boundary, and therefore, to an increased resistance to plastic deformation being responsible for the delay of the shear localization process.

Figure 10

Sample S₁, length scale-dependent behavior: spatial distribution of free volume ξ at an applied compression strain |E₂₂|=8% for different

ratios of internal length l and external specimen size l₀.

In summary, a larger

ratio first leads to a change from catastrophic shear localization, which is usually observed in bulk metallic glasses, to stable shear band formation. A further increase causes a delay and finally even a suppression of the shear band formation.

4.3 Rate-dependent behavior

In this section, the rate dependence of the model is analyzed. As known from the phenomena related to microstructure evolution, the strain rate has a strong influence on the material behavior and may lead to change and also the suppression of the microstructure and its development [61, 74, 75].

Figure 11 displays the average stress-strain behavior of the model for four different strain rate ratios

of external loading rate

to internal material rate

As expected, no effect is observed in the elastic region. As usual, with increasing ratio

the stress level increases and yielding occurs at a larger strain. Therefore, the highest ratio

always shows the highest stress. The general model behavior is the same for all rates; however, small differences arise: after yielding, the negative slope of the stress plateau increases with increasing

Figure 11

Sample S₁, rate-dependent behavior: average stress component |σ₂₂| over applied strain |E₂₂| under compression loading for different

ratios of external loading rate

to internal material rate

The stress drop is steeper and longer until the saturation value is reached for larger rate ratios. With ongoing deformation, the stress difference between the results is staying constantly for the different rates.

The spatial distributions of ξ, c, and γ are very similar for the different rates. For lower strain rates, a slightly more diffuse distribution, which gets sharper with increasing rate, is observed. The results show that our model for metallic glasses behaves similar to most crystalline materials with respect to its rate dependence. An increase in the loading rate generally leads to a stress increase.

4.4 Tension-compression asymmetry

In the following part, the tension-compression asymmetry of our model is investigated. The asymmetry of metallic glasses under tensile and compressive loading is one of its distinguishing attributes compared with other materials such as metals. Metallic glasses start to yield earlier under tension than under compression (see, e.g., [17, 76–78]). To investigate this behavior, Figure 12 shows the average stress-strain behavior under tension and compression. As expected, an earlier yielding is observed for tension. With ongoing deformation, the tensile behavior is similar to the compressive behavior; however, the knee region disappears (similar as observed for increasing

).

Figure 12

Sample S₁, tension-compression asymmetry: average stress component |σ₂₂| over applied strain |E₂₂| under tensile and compressive loading.

Figure 13 shows the spatial distribution of the free volume ξ under compressive and tensile loading. As observed, under tensile loading, the sample shows a necking behavior instead of the formation of a bulge structure. In addition, the distinct shear band is more pronounced under compression, whereas the side shear band is more pronounced under tension.

Figure 13

Sample S₁, tension-compression asymmetry: spatial distribution of free volume ξ at an applied strain |E₂₂|=8% under compressive (left) and tensile (right) loading. The shear band angle under compression is Θ=41° and under tension Θ=46°.

The shear band develops at a different angle under tension compared with compression. In case of compression, a shear band angle of Θ=41° is observed in the simulations. This is in excellent agreement with experimental findings of Wright et al. [15], who observed a compressive fracture angle of Θ=42° as well as with the theoretical study of Schuh and Lund [17], who obtained Θ=41.5° for the compressive shear angle. In the case of tension, the shear band is formed at an angle Θ=46°. In the experiments, a fracture angle Θ between 50° and 56° is reported under tension (e.g., [15, 78, 79]). However, these studies were performed for bulk metallic glasses and size effects were not explicitly investigated. These studies also report that the bulk metallic glass does not yield but rather shows a direct brittle fracture behavior under tension. Under compression, a catastrophic failure behavior is reported, which changes with decreasing size [26] (cf. Section 4.2). One explanation for the different behavior under tension and compression was given by Jiang et al. [80], who proposed that the mechanism of direct breakage of local atomic clusters involves a tension transformation zone next to the shear transformation zone under tension. If the material does not have enough time to fully flow and therefore relax, the shear transformation zones is restrained and fracture occurs via the tension transformation zone. For more details, the reader is referred to [16, 80].

In summary, the model is able to predict the tension-compression asymmetry accurate for specimens of small size. The shear angle for compressive loading fits quite well with experimental observations. However, for tensile loading, a small deviation is obtained with respect to experiments of large-sized bulk metallic glass samples. This might be due to size effects.

4.5 Investigation of the influence of initial defects on the example of sample S₂

In this section, the model behavior for sample S₂ is studied. This sample represents a square loaded in tension with periodic deformation and free-volume ξ boundary conditions. Based on this numerical example, the influence on the position of the initial defect, represented by a lower initial cohesion c₀, is studied. Figure 14 shows the spatial distribution of the free volume ξ at three applied tensile strain states for a defect near the center of the structure. In the beginning, a similar distribution as for sample S₁ is obtained (see Figure 4), but with increasing strain, differences occur. These are related to the different structure and periodic boundary conditions. Due to the periodic free-volume boundary conditions, parts of three shear bands are observed within the sample. Sample S₂ leads to the same response under tensile as well as compressive loading due to the periodic boundary conditions, which, e.g., suppresses necking phenomena. In this example, the shear band forms from the bottom left corner to the top right corner. The direction is defined by the position of the defect. Depending on this position, different distinct shear bands are observed. The cohesion c and the plastic strain γ are showing a similar behavior as the free volume ξ.

Figure 14

Sample S₂, variation of initial defect: spatial distribution of free volume ξ at different applied strains E₁₁.

Initial defect near the center (see Figure 16).

Next, the spatial distribution of the free volume ξ at three applied tensile strain states for a defect at the center of the right edge is shown in Figure 15 (see also Figure 16). In this case, the shear band proceeds from the bottom right to the top left due to the fact that the initial defect is located above the horizontal center axis.

Figure 15

Sample S₂, variation of initial defect: spatial distribution of free volume ξ at different applied strains E₁₁.

Initial defect at the right edge slightly above horizontal center (see Figure 16).

Figure 16

Sample S₂, variation of initial defect: spatial distribution of free volume ξ at an applied tensile strain E₁₁=8% for different places of the initial defect.

In each simulation, only one element (shown in black of the mesh in the middle) had a prescribed lower initial cohesion of 0.99c₀ as indicated.

To summarize the influence of the position of the defect, Figure 16 displays the results of five simulations with different positions of initial defect as indicated in the figure. Differences are observed with respect to the distributions and the directions of the shear band. The average stress-strain behavior is the same for all simulations. Although different defect positions induce different nucleation position of the shear bands, the RVE predicts the same overall material behavior.

Finally, two limiting cases are investigated: (i) randomly distributed initial cohesion {0.99–1}c₀ and (ii) a perfect structure without initial defect. Figure 17 depicts the free-volume distribution ξ. For case (i), a much sharper shear band is obtained compared with the previous simulations, which is due to a different initial kinetic state. The cohesion c has a kinetic character and does not change the internal energy state. The heterogeneous initial kinetic state leads to an imbalanced atomistic structure. The atoms aim relaxation and at reaching the lowest energy state. As a consequence, a much sharper shear band develops faster. This finding is visible in the average stress-strain curve in Figure 18. The larger number of defects leads to an acceleration of the shear localization process.

Figure 17

Sample S₂, variation of initial defect: spatial distribution of free volume ξ at an applied tensile strain E₁₁=8% for sample with randomly distributed defects (left) and sample without defects (right).

Figure 18

Sample S₂, variation of initial defect: average stress component |σ₁₁| over applied strain |E₁₁| under tensile loading for different conditions for the initial defects within the sample.

Solid line, one element with defect placed near the center; dashed line, results with randomly distributed defects; dot-dashed line, no defect placed in the structure.

As a last example, case (ii) is investigated. The free-volume ξ distribution at 8% strain is shown in Figure 17 (right). Although no defects are existing in this sample to trigger the shear band formation, shear bands that start from the center of each edge are observed. However, the shear bands form slower than in case of a sample with defects. Therefore, a pronounced shear band cannot yet be observed at E₁₁=8%. Another consequence is a slightly higher stress level once a shear band pronounces in the samples with defects (see Figure 18).

The change of the defect’s position to trigger the shear band formation does not affect the macroscopic behavior and the shear band structure. However, it affects the shear bands direction but not its orientation. Increasing the number of defects leads to an imbalance of the atomic structure and, consequently, to an acceleration of the shear band formation. Without internal defects, shear localization is observed, however, slightly delayed.

5 Summary

In this work, a nonlocal thermodynamically consistent model formulation is presented for metallic glasses. The numerical implementation is carried out with the help of the finite element method. Based on several numerical experiments, the model behavior was investigated. The presented model formulation is able to model the stable growth of shear localization in submicron samples. A length scale-dependent behavior is observed. The size effect leads to a delay and, finally, even a suppression of shear localization. An increased loading rate leads to higher stresses in the metallic glasses. The model also maps the tension-compression asymmetry. The shear band angle under compression is predicted as Θ=41° and as Θ=46° under tension. Further numerical experiments were performed on the basis of a periodic sample with multiple shear bands. The influence of the defects to trigger the shear band formation was investigated. The position of the initial defect(s) does not affect the macroscopic behavior – only the direction of the shear band evolution. It is shown that the inclusion of several defects leads to an imbalance in the atomic structure resulting in an acceleration of the shear band formation.