A nonlocal method for modeling interfaces: Numerical simulation of decohesion and sliding at grain boundaries

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 362 (2020) 112836 www.elsevier.com/locate/cma

A nonlocal method for modeling interfaces: Numerical simulation of decohesion and sliding at grain boundaries Shahed Rezaeia ,∗, Jaber Rezaei Mianroodib,c , Kavan Khaledia , Stefanie Reesea a

Institute of Applied Mechanics, RWTH Aachen University, D-52074 Aachen, Germany b Material Mechanics, RWTH Aachen University, 52074 Aachen, Germany c Microstructure Physics and Alloy Design, Max-Planck-Institut für Eisenforschung, Düsseldorf, Germany Received 10 April 2019; received in revised form 18 December 2019; accepted 7 January 2020 Available online xxxx

Abstract Understanding and modeling the interface behavior is an important task for predicting materials response in various applications. To formulate the behavior of an arbitrary interface, one needs to construct the relation between acting tractions and displacement jumps at the interface. In addition to capturing the correct physics of the interface, the so-called traction– separation relation must also be thermodynamically consistent and satisfy the basic balance laws. Apart from many attempts in the literature to address these issues, a new and simple method to capture the complex mechanical behavior at an arbitrary interface is proposed. The new formulation is based on introducing a new quantity called “traction density". As a result, the traction–separation relation for any arbitrary interface is automatically computed by integrating the traction density over the interface. The traction density can be formulated based on understandings and observations from lower scales. As will be shown, the mathematical representation of the traction density is relatively simple and therefore its consistency can be verified easily. When it comes to the grain boundary (GB) behavior, the proposed methodology is able to represent not only intergranular fracture but also grain boundary sliding. For calibration and verification of the model, molecular dynamics (MD) simulations for aluminum Σ 5 GB are utilized. Interestingly, the calculations from current MD simulations show size-dependent behavior for the GB. By introducing a healing parameter in the new interface model, it is now possible to explain and predict possible GB size-dependent behavior. c 2020 Elsevier B.V. All rights reserved. ⃝ Keywords: Cohesive zone element; Molecular dynamics; Interface plasticity; Intergranular fracture; Healing; Size-dependent TSL

1. Introduction Understanding interface mechanics plays a major role in determining and predicting the overall material behavior in many applications. The grain boundary is an interesting and also important example of an interface, which influences the fracture resistance of many types of materials [1,2]. The contribution of the GB is more pronounced when it comes to noncrystalline metals and ceramics [3,4]. Earlier investigations by means of experiments or atomistic simulations show that the GB behavior is in general anisotropic. In other words, by loading the GB in ∗ Corresponding author.

E-mail address: [email protected] (S. Rezaei). https://doi.org/10.1016/j.cma.2020.112836 c 2020 Elsevier B.V. All rights reserved. 0045-7825/⃝

2

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836

Fig. 1. Traction for different interfaces during loading and unloading.

different directions, various mechanical behaviors might be activated individually or simultaneously. As an example, studies by Rezaei et al. [5] on aluminum GB show intergranular fracture being more pronounced during loading normal to the GB (mode I) than for mode II or III. On the other hand, GB sliding is mainly reported for shear loading (mode II or III). The anisotropic behavior of the interface is also reported in other applications such as delamination of composites [6]. It is a well-established approach to utilize an interface model to include and capture the aforementioned points. The interface model is usually implemented in the finite element (FE) context by means of an initially zero thickness element also known as cohesive zone (CZ) element [7]. Later on, by means of enriching all the interfaces via the interface elements, one is able to capture the interface behavior and predict its influence on the overall results [8,9]. Despite many developments for the interface modeling [10–13], one remaining and challenging aspect is to physically relate the traction vector to the displacement jumps or discontinuities [14]. This includes proper coupling of different fracture modes [15] and also plasticity and damage at the interface [16]. Additionally, these models have to be thermodynamically consistent (see [17–19]). Furthermore, the interface model should satisfy the basic balance laws such as the balance of angular momentum. This holds in particular, when it comes to finite displacements and large deformation (see [20] and [21]). Many attempts have been made to overcome the mentioned problems and to improve the formulation of the interface elements over the past decades. Bosch et al. [22] investigated the delamination in polymer coated steel by proposing an extended cohesive zone model for large displacements which allows for a mode-dependent behavior in the large deformation case (see also [23]). Vossen et al. [24] showed that for many existing CZ models, considering anisotropic behavior would result into violation of the balance of angular momentum at the cohesive element level. In other words, to satisfy the balance of angular momentum, the traction vector should be aligned with the gap vector. The error introduced by this assumption might be negligible, if the ratio between the critical opening length and the length of the process zone is small. In order to capture anisotropic behavior of the interface at large interface opening in a thermodynamically consistent way, Mosler and Scheider [25] showed the necessity of considering additional stress tensors (deformational surface shear). A similar approach is also taken by Ottosen et al. [26]. An efficient finite element implementation of this framework is discussed by Heitbreder et al. [27]. Distribution of tractions during loading an arbitrary interface is shown in Fig. 1. On the left-hand side, an isotropic relation for the traction–separation law is shown. The former assumption is often utilized depending on the application and required accuracy. On the other hand, on the right-hand side of Fig. 1, other possible complexities at a general interface are shown. This includes the anisotropic behavior of the traction vector (note that tractions are no longer colinear with the gap vector) and also the additional remaining gap upon unloading. The main questions regarding interface modeling can be summarized as follow. How can we have an anisotropic relation for the traction–separation law and still satisfy the balance of angular momentum? How can the decohesion process together with sliding be modeled accurately? What would be the required internal variables and evolution laws? Moreover, models should satisfy the basic balance laws and be thermodynamically consistent. It is also a big advantage if the developed models are easy to implement and do not require so many input parameters. Although these important questions have been addressed partially in previous works [14,26,28,29], there are still some missing points. For example, what is the underlying physics behind the anisotropic behavior? Or apart from

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836

3

satisfying the balances of linear momentum and energy, how different mechanisms such as plasticity and damage should be coupled? We strive to answer these questions by understanding the physics of each particular interface at smaller scales and use this knowledge at the upper scale. This is similar to the idea upon which the FE2 method is based, where the macroscopic stress–strain behavior is obtained from computations at the micro-scale. As will be shown later, considering the nonlocal nature of the traction vector at the interface can help us to solve the above mentioned problems. The idea of considering nonlocal interaction has been addressed before by many authors [30–33]. Recently a phase-field approach is also proposed for interface failure in which an internal length scale is considered in the formulation [34] (see also [35,36] for an extension to gradient-based nonlocal damage). From the discretization point of view, one of the drawbacks of using the classical CZ elements is the fact that the meshes (positions of the nodes) at the opposite sides of the interface must initially match each other. Few extensions have been introduced in the literature using the discontinuous Galerkin method [37] and the isogeometric CZ modeling [38] to solve this problem. An interesting formulation which overcomes the problem of having nonmatching meshes is introduced by Paggi and Wriggers [29]. The introduced node-to-segment CZ element seems to be capable of (additionally) solving the problem with satisfying the balance of angular momentum. Utilizing such a discretization approach seems to be promising in introducing interface elements for meshing interfaces in complicated structures. Determining the interface parameters is another challenging aspect, especially when it comes to simulations at small scales [39]. Utilizing atomistic simulations was shown to be a good instrument to obtain parameters for continuum interface models (see [40–44] to name only a few). Transferring information from small scale to large scale is also an important issue [45]. Besides the parameter identification, atomistic simulations also provide a deeper insight into different active mechanisms at the interface [46,47]. Elzas et al. [48] investigated the GB interface under shear loading. Their studies show that the behavior of the GB in shear direction (shear traction versus shear separation) is very similar to elastoplastic continuum models with hardening (see also [5]). Furthermore, through the experiments on finite size grains, a size effect was observed for GB sliding [49]. To the authors’ best knowledge, the latter observation was never verified and investigated numerically. The sliding of grains along the GB is a result of two simultaneous mechanisms, i.e. (1) the loss of the atomic bond between the two neighboring atoms and (2) the formation of new atomic bonds in the direction of sliding. The latter point is related to healing phenomena which calls for more attention into the modeling of grain boundaries. Atomic simulations also show that healing happens, whenever the distance between the two interfaces becomes smaller than a certain value [5]. The healing effects or re-bonding of atomic structure have a major influence on the overall material behavior. Joining similar or different types of metals (e.g. by cold welding) is also related to the recovery of atomic bonds at small scales [50–53]. From the modeling point of view, taking into account the healing effects is an interesting and challenging task on its own (see [54]). Healing mechanisms were included into cohesive zone formulations by only a few authors [55,56]. Therefore, the physical and consistent introduction of the healing parameter needs more attention and more experimental evidence. The main idea here is to obtain the traction at the interface by means of the information from lower scales. The interface under focus in this work is the grain boundary. Nevertheless, it should be possible to extend the proposed methodology to other types of interfaces. The structure of the paper is as follows. In Section 2 the main problem and motivation is described. In Section 3 the new proposed methodology and the connection between traction and traction density is discussed in detail. Section 4 is dedicated to the formulation behind traction density for different applications such as grain boundaries. Moreover, the finite element formulation and discretization of the new interface element is explained in Section 5. The results are finally reported in Section 6, where for the calibration of the model and also further verification a comparison with calculations from MD simulations is carried out. Finally, in Section 7 a discussion is provided for the readers about the potential of the current methodology and its comparison with other interface models in the literature. The future works and challenges are addressed later on as well. 2. Problem description The configuration of a general body B and its boundary conditions are shown in Fig. 2. Without losing generality, the focus is here on quasi-static problems (no inertia forces). The prescribed external traction and displacement vectors are denoted by t p and u p , respectively. The governing local form of the equilibrium equations in the current

4

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836

Fig. 2. Prescribed displacements and tractions on a general body (left). Formation of a crack and cohesive tractions (right).

configuration is summarized in Eqs. (1)–(3). In the present modeling approach, for simplicity, it is assumed that the body B can be described by these classical local equations (i.e. it does not include any nonlocal effects). ∇ ·σ +b=0

on

u = up σ · n = tp

B = B+ ∪ B−,

(1)

on Γu ,

(2)

on Γt ,

(3)

Here, σ is the Cauchy stress tensor, b represents the vector of body forces, u is the displacement vector and n is the normal vector to the surface Γt . Moreover, although not pursued in this work, it might be more suitable to use nonlocal models within the body B as well. As an example, one can mention the micromorphic approach developed by Forest (2009) [57] in which the nonlocal interactions are taken into account by introducing micromorphic variables (see also [58]). Note that in the current work, the continuum model for the body B, is local, and the nonlocal effects are only limited to the predefined interface, as it is discussed in the next section. The interface Γc defined in the body B represents the location of a possible displacement discontinuity. After applying the boundary conditions, the original body B is divided into two bodies B + and B − . The interface Γc is divided into two parts as well, and tractions t + and t − should be considered along with the two new interfaces. The associated normal vectors are n+ and n− , respectively. The effect of the cracked region is only treated as a traction boundary condition on the classical body, B, through Eqs. (4)–(5). σ · n+ = t +

on Γc+ ,

(4)

σ · n− = t −

on Γc− .

(5)

The idea behind the different interface models is to formulate a relation between the traction vector t and the displacement jump vector JxK at the interface (surface Γc ). On the other hand, finding a general relation is a challenging task especially when it comes to large displacement jumps at the interface. The idea of the present work is to introduce an alternative method to alleviate these problems. 3. Concept of traction and traction density The idea behind current work is to obtain the traction vector t by integrating the values of the traction density vector T . This procedure and formulating the weak form is explained in the following. In the left side of Fig. 3, X + represents the reference position of an arbitrary point i belonging to body B + at the interface. The position vector X − represents the position of the point i projected on the lower interface belonging to body B − . In the undeformed state, these two vectors are equal as shown in Fig. 3. The vector x + (x − ) represents the current position of the point i on the upper (lower) surface Γ + (Γ − ). The so-called gap vector related to point i is defined by g i = x i+ − x i− . Further, we define for each point x i− on the lower surface Γ − the so-called local gap vector G i (x − ) which represents the vector from an arbitrary point on Γ − to the point i. The distribution of the traction density vector related to point i is represented by T i (x − ) (see green arrows and region on the right hand side of Fig. 3).

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836

5

Fig. 3. Reference and current state of the interface. The assumption of collinearity of T and G is made here. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

In general, the traction density vector T is a function of the current local gap between the two arbitrary points and a set of internal variables v, i.e. T i = T i (G, v).

(6)

Note that for further simplification, influences of other parameters such as temperature or loading rates are excluded in this work. Furthermore, in the current work, the focus is on mechanics of the GB and as it will be shown in the following section, the only internal variable will be the scalar damage variable D. Considering other applications such as fiber bridging in delamination of composites, other variables may be added in the traction density formulation. Finally, the traction vector at point i can be obtained by summing up the traction density on the lower surface ∫ ti = T dΓ − . (7) Γ−

What is shown via Eq. (7) is the key idea behind the current work. Here, the traction vector t is not introduced as a direct function of the gap vector g like the conventional CZ models. Instead, it is obtained based on the traction density vector T which depends on the local gap vector G. By integrating the traction density, the traction vector g automatically includes some nonlocal information. As a result, the global traction vector t is no longer an isotropic function of the global gap vector g. As it is shown in Fig. 4, the traction vector at the interface is determined by summing up the traction density vector in the vicinity of each point. The area (segment) on Γ − which has an influence on Γ + (or other way around), depends on the nature of the problem. Here, it is assumed that the two surfaces (Γ + and Γ − ) attract each other by means of atomistic forces, which explains why the traction density vanishes eventually as the distance between the interfaces increases. Therefore, for each arbitrary point i shown in Fig. 4, the dashed lines represent the area on the lower interface Γ − (or Γ + ) which influences (attracts) the point i. The traction vector t + and t − in Eqs. (4) and (5) can be written in terms of the traction density vector T + and − T , respectively: ∫ + + t (x ) = T − (x + , x − , v) dΓ − , (8) − Γ ∫ t − (x − ) = T + (x + , x − , v) dΓ + . (9) Γ+

Balance of linear and angular momentum must hold within the interface. Therefore, one possible choice would be to assume that the traction density vector from opposite sides of the interface are equal in magnitude and also co-linear with the local gap vector G = x + − x − (see also Appendix B) T − = −T + = T , G ∥ T ⇒ G × T = 0.

(10) (11)

6

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836

Fig. 4. Relation between traction and traction density vectors.

Since the traction vector t + (t − ) is computed by integrating the vector T + (T − ) according to Eq. (8) (Eq. (9)), it also satisfies the balance of angular momentum. Having developed the relation for the traction vector t and taking into account equations (1)–(5), the weak form is given by ∫ ∫ ∫ ∫ − σ : ∇δu d V + t + · δu+ dΓ + + t − · δu− dΓ − + t p · δu dΓt = 0. (12) + − B Γt Γ  Γ  Icz

The contribution of the cohesive virtual work Icz is further elaborated using Eqs. (8)–(10) ∫ ∫ ∫ ∫ T − (x + , x − , v) · δu+ dΓ − dΓ + + T + (x + , x − , v) · δu− dΓ + dΓ − Icz = + Γ− − Γ+ Γ Γ ∫ ∫ ∫ ∫ = T · δu+ − T · δu− dΓ − dΓ + = T · δG dΓ − dΓ + . Γ+

Γ−

Γ+

(13)

Γ−

4. Formulation of traction density 4.1. Study of the grain boundary The idea is to look at the interface at a smaller scale down to the atom level. The interface at different scales is shown in Fig. 5. On the left side of Fig. 5, where the sizes are among a few Angstrom, the separation of two atoms is shown. The force F acting between two atoms strongly depends on the distance r (apart from other parameters such as possible electrical charge). This relation is known to be completely reversible. On the right side of Fig. 5, where the dimensions are about a few micrometers, the traction vector t is shown in terms of its components in normal (tn ) and shear direction (ts ). The interface traction, in general, is of anisotropic nature and shows different behavior in different directions. The latter point is also clear from different measured values for the fracture toughness or maximum traction after mode I and mode II loading [5,59]. An intermediate scale is introduced between the atomic scale and microscale. The physical dimensions of this scale can be considered to be a few nanometers (see middle part of Fig. 5 shown in the light green box). Considering a block of atoms from the upper and lower grain, the traction density vector T can be defined as a quantity which describes the traction between the two interfaces per unit of area. The introduced traction density function has two important features: • Due to the fact that the blocks of atoms belonging to two neighboring grains are small enough (and ignoring other physical aspects such as electrical charges and possible chemical effects), the traction density vector T is assumed to be co-linear to the local gap vector G (T ∥ G). • The blocks of atoms are large enough so that they include structural features such as damage or reduction in the initial interface stiffness. In other words, the relation between ∥T ∥ and ∥G∥ is not completely reversible and includes some history of the system here expressed by means of the damage variable D.

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836

7

Fig. 5. Separation of two neighboring atoms (left). Separation of colony of atoms from opposite side of the grain boundary (middle). Traction vector acting on the grain boundary at microscale (right). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The proposed relation for the traction density T versus the distance between two blocks of atoms G is described via Eq. (14) and is shown in the left part of Fig. 6. Within the current work, the parameters to describe the traction density are as follows: the maximum traction density (T0 ), the distance at which the maximum traction density is reached (∆0 ), the distance at which the traction density vanishes (∆ f ) and the healing distance parameters (Hl ). The traction density vector is written in the form T = (1 − D) K 0 G,

(14)

where K 0 = T0 /∆0 and the scalar damage variable D sets in, as soon as the traction density reaches T0 or the opening passes ∆0 (see also [13] and [12]). Therefore, the reduced√stiffness is defined as K D = (1 − D)K 0 . The distance between the two interfaces is obtained by ∥G∥ = ∆ = ⟨G n ⟩2 + G 2s and the damage variable is going to be determined by the following simple steps. Considering the healing effects, the damage variable D may decrease, if the distance ∆ becomes smaller than the healing length Hl . Therefore, the variable ∆m is introduced: { max(∆) if ∆ ≥ Hl ∆m := (15) ∆ if ∆ < Hl Here ∆m denotes the maximum distance achieved during the loading when the distance ∆ is greater than the healing distance Hl . The value of ∆m resets as soon as healing happens during the unloading of the system (i.e. ∆ < Hl ).

8

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836

Fig. 6. The behavior of traction density (left) and the damage variable (right) as a function of the gap for the grain boundary.

Finally, the damage variable D is ⎧ if ⎪ ⎨0 ∆f ∆m −∆0 if D = ∆ f −∆0 ∆m ⎪ ⎩ 1 if

defined as follows (see also right hand side of Fig. 6): ∆m ≤ ∆0 ∆0 < ∆m < ∆ f

(16)

∆ f ≤ ∆m

The healing length Hl is usually selected to be very small (0 < Hl < ∆0 ). The healing effect increases the damaged stiffness K D during the unloading. The formulation introduced at this point is motivated from the MD simulations. Physically, the traction density for the GB can be interpreted as a long-range bond between material points, analogous to the interaction force between atoms. Indeed, even initially, most of the interaction forces (and traction density vectors) are non-zero. Therefore, the system includes eigenstresses before loading. This picture is in complete agreement with atomistic systems, where near defects (such as cracks) the lattice structure is disrupted and distorted, creating non-zero but balanced interaction forces between atoms. Note that modeling macroscopic systems with fully resolved atomistic details are beyond the current computation power. The current model of traction–separation law for material point interactions is non-conservative, through Eqs. (14), (15), (16). Thus, the model is able to homogenize the atomistic details leading to irreversible damage at the continuum scale, without explicitly tracking the atomistic details of fracture surface or environmental contamination. At this point, more advanced phenomenological models may be utilized to address the healing effects (e.g. see [54,60]). Doing so, additional parameters may be introduced to control the transition between damage and healing in a regularized way. Furthermore, the relation in Eq. (16) is based on the bi-linear traction density relation and can be further extended depending on the application. Consistency of the formulation can be also elaborate considering the proposed formulation for the traction density. During the loading, the dissipation is always positive through Eqs. (14), (15), (16) (see also [8]). Further derivation should be covered in future works. It is worth to mention that the proposed model is checked for various cases against the atomistic results (see Section 6.3). 4.2. Traction density in further applications The introduced method is applicable to interfaces at a completely different scale. An example is inter-laminar failure and damage progression [61,62] as shown in Fig. 7. Here, the upper surface is connected to the lower surface by means of a complex fibers distribution. Again, the influenced area and the relation for the traction density can be determined based on the fiber distribution and the mechanical behavior of each individual fiber [63]. The traction density here is simply co-linear with the fiber direction or the local gap opening G. On the other hand, by integrating the traction density over the interface, anisotropic traction will come out naturally. See also studies on intra-laminar fiber bridging in composites by Sørensen and Jacobsen [64]. 5. Finite element implementation The discretization of the interface and the contribution to the residual vector are discussed. The formulation is presented for the 2D case and the extension to the 3D case is straightforward.

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836

9

Fig. 7. Representation of the fiber bridging via traction density in delamination of composites.

Fig. 8. Discretization of the interface using node to segment elements.

In the left side of Fig. 8, the finite element mesh for the body B in the reference configuration is shown. As an example, one interface element is denoted by its nodes 1, 2 and 3. The interface element’s thickness is initially zero. The first two nodes lie on the lower surface, while the third node lies on the upper surface. Point 3′ only denotes the initial position of the third node once the interface is closed and does not represent any additional node. According to the node numbering described in Fig. 8, the reference position vector X vT , the current position vector x vT and the displacement vector uvT of the interface element nodes can be written in vectorial form [ y y y] X vT = X 1x X 1 X 2x X 2 X 3x X 3 , [ ] y y y uvT = u 1x u 1 u 2x u 2 u 3x u 3 , x vT = X vT + uvT .

(17)

The gap vector between an arbitrary point on the lower surface and node number 3 is given by G(ξ ) = R Jx(ξ )K = R B(ξ ) x v .

(18)

As shown in the middle part of Fig. 8, the coordinate ξ in the parent element varies between −1 and +1. Furthermore, R is the rotation matrix which is defined based on the rotation of the lower interface compared to the reference coordinates X and Y (see angle θ in the right-hand side of Fig. 8): [ ] cos (θ) sin (θ ) R= . (19) −sin (θ ) cos (θ ) For the current formulation, introducing the rotation matrix is not needed. This also can be considered to be an advantage of the formulation. Nevertheless, since for practical application one has to take care of the contact situation to avoid normal penetration of the surfaces, this part is added. Here, the negative normal gap is penalized via a

10 S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836

simple penalty approach (see [8] for more details). In addition, by means of the matrix [ ] −N1 0 −N2 0 1 0 B= , 0 −N1 0 −N2 0 1

(20)

the current gap G is interpolated through the interface from the nodal values x v , where linear shape functions N1 = 0.5 (1 − ξ ) and N2 = 0.5 (1 + ξ ) are utilized for the interpolation. The overall gap vector for each interface element is also defined as g = G(ξ P ),

(21)

which, as shown in Fig. 8, measures the distance between node 3 and its initial position on the lower segment 3′ . Here, ξ P represents the coordinate of node 3′ in the parent element, i.e. x 3′ = N1 (ξ P )x 1 + N2 (ξ P )x 2 . Depending on the nature of the problem the traction density values at the interface can vary in a completely nonlinear manner (see the middle part of Fig. 8). Therefore, as discussed in Appendix A, for the current formulation it is important to consider enough integration points on the segment 1–2 in order to be able to capture the physics of the interface accurately. Furthermore, the traction vector at node 3 (see Fig. 8) can be obtained by summing up the contribution of traction density from different elements (see also Eq. (7)) ) nI P ( ∑∑ L− el wi . (22) t3 = Ti 2 e e i=1 The interface contribution in the weak form can then be evaluated as ) ( ∫ ∫ ∫ ∂G T − + T T dΓ − L + Icz = δG · T dΓ dΓ = δuv el ∂ uv Γ+ Γ− Γ− nI P ∑ ( T T ) L− el = δuvT B i R T i wi L + , el 2 i=1

(23)

+ where L − el is the length of the interface element or the segment 1–2 and L el is the sum of half the distance between node 3 and its neighboring nodes on the upper surface (see Fig. 8). Considering Eqs. (18) and (21), the notation B P = B(ξ P ) is used in Eq. (23). Furthermore, wi is the weighting for the ith integration point and n I P represents the number of integration points. According to Fig. 6 and Eq. (14), the following relation holds for the traction density:

T i = (1 − Di ) K 0 G i = (1 − Di ) K 0 R B i x v .

(24)

Finally, the linearization of Icz (see Eq. (23)) with respect to the nodal displacement vector uv is carried out (see also [65]): (n ) IP ∑ L− el T T ∆Icz = δuv (R B P ) C i R B i wi L+ (25) el , 2 i=1 The derivative of the traction density vector with respect to the gap vector is obtained by dT ∂T ∂T ∂ D ∂T ∂ D C= = + = (1 − D) K 0 I + . (26) dG ∂G ∂D ∂G ∂D ∂G Depending on the healing parameter Hl , the area from the lower surface which has an influence on a sample point from the upper surface may constantly change. Therefore, for the mesh generation procedure, one should take into account additional interface elements which connect the nodes from one surface to segments far away from the lower surface (see Fig. 9). Utilizing the node-to-segment type of discretization for the interface, one should make sure to use enough number of elements to assure symmetry of results and reduce discretization error as much as possible. Certainly, other element types (e.g. conventional four-node elements) can also be applied to discretize the interface with the current formulation. At this point, one may also notice similarities between the proposed methodology and the basics of peridynamics which is an alternative method for continuum mechanics. According to Silling et al. [33,66,67], in peridynamics, the interaction of different nodes is defined in terms of a bonding which acts in a certain horizon. See also [68,69]

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836 11

Fig. 9. The active region at the interface during an arbitrary interface opening. Table 1 Properties of the introduced interface model for parameter studies. Maximum traction density T0 Opening at the maximum traction ∆0 Critical interface opening ∆ f Healing length parameter Hl

Unit

Value

[GPa/µm2 ] [µm] [µm] [µm]

10.0 5 × 10−2 2 × 10−1 2 × 10−2

and [70] for application of peridynamics in the simulation of fracture in polycrystalline materials. Another interesting extension could be considering the contribution of atoms below the surface. At the moment their influence is lumped into the formulation of the traction density vector. For further development, one may consider different healing distances depend on the separation mode. 6. Results 6.1. Traction–separation relation In the following, the introduced formulation is applied at the interface between two (rigid) grains as shown in Fig. 10. For mode I loading, the grains are separated in Y direction (or direction of n) and for mode II they are sheared against each other in X direction (or direction of s). The interface size is denoted by L x . The computations are based on the plane strain assumption and the out-of-plane thickness of the model is L z = 1 µm. Some possible finite element meshes are shown in the lower part of Fig. 10 to make more clear how the nodes from the opposite interface are connected to each other. For the current study, the upper and the lower interface are discretized by 100 elements. The interface parameters are given in Table 1 unless mentioned otherwise. The grains are considered to behave very stiffly in comparison with the interface. The displacement of the upper edge of the grain 1 increases proportional to the pseudo-time step (no rate dependency is considered for now). The values for the normal and shear tractions reported in the following are obtained by measuring the reaction forces Fx and Fy at the top edge divided by the length L x = 1 µm. The influence of the interface maximum traction density is studied in Fig. 11. For a better comparison, the system is also unloaded at some point. Note that the dashed lines in Fig. 11 represent the original loading path and not the reloading path. There are two important features about the current formulation. First, the unloading path for mode I is similar to elastic damage models in which the overall normal traction goes back to zero without any permanent normal gap. Moreover, due to healing effects, the normal traction increases as soon as the interface gap becomes smaller than the healing length (gn < Hh = 0.02 µm). Secondly, for mode II loading, a permanent shear gap is predicted as the system is unloaded. Please note that no plasticity formulation or additional internal variables are added to the formulation. The only internal variable is the scalar damage variable D which represents degradation between atomic bonds. The healing distance Hl (a constant parameter) plays an important role in capturing the interface plasticity. Simply by considering the influence of damage and healing of the traction density T , many important features are captured automatically. The influence of the parameter ∆0 is studied in Fig. 12. As expected, for mode I loading, the initial stiffness of the interface reduces by increasing ∆0 . On the other hand, the shear tractions are less affected by this parameter. The influence of the critical interface opening ∆ f is studied in Fig. 13. Considering mode I and mode II loading, increasing ∆ f results into higher maximum traction at the interface as well as maximum final gap before the full failure.

12 S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836

Fig. 10. Introduced interface model under different loading cases and its discretization.

Fig. 11. Influence of the parameter T0 while keeping other interface parameters constant (see Table 1).

The influence of the healing length parameter Hl is studied in Fig. 14. By increasing Hl the amount of gain in the traction through the unloading for mode I increases. For mode II loading, the value of the shear traction only vanishes as the shear gap becomes larger than the interface length and critical opening, i.e. gs > L x + ∆ f . If the healing distance is reduced to a very low value (Hl = 0.001 µm), the influence of the healing vanishes. By decreasing the healing distance Hl , less plastic gap (remaining shear gap) is observed (Compare the unloading slope in Fig. 14(b)).

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836 13

Fig. 12. Influence of the parameter ∆0 while keeping other interface parameters constant (see Table 1).

Fig. 13. Influence of the parameter ∆ f while keeping other interface parameters constant (see Table 1).

The results for mixed-mode behavior of the interface model are shown in Fig. 15 for different separation angles φ (see Fig. 10). 6.2. Influence of the interface length In this section, size-dependent traction–separation relation is studied. For this purpose, the interface length L x is changed and the obtained normal and shear tractions are reported in Fig. 16. In mode I loading, the interface length barely changes the results for the normal traction. On the other hand, mode II results are strongly dependent on the interface length. This is due to healing effects and formation of new bonding as the two interfaces slide on each other. It should be noted that in the case of Hl = 0.0 (when the healing effects are ignored), the interface behavior is almost independent of the interface size L x . According to Fig. 17, the overall traction–separation relation in both loading modes converges to certain values when L x is large enough. Note that the behavior of the material points near the free surfaces at x = 0 and x = L x is influenced by the boundary conditions. With increasing length of L x , the boundary conditions influence an increasingly smaller number of material points (see Fig. 17).

14 S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836

Fig. 14. Influence of the parameter Hl while keeping other interface parameters constant (see Table 1).

Fig. 15. Obtained tractions for different values of angle φ while keeping other parameters constant. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

6.3. Comparison with MD simulation The main motivation for the following section is to compare the obtained traction–separation relation from the model with the results of the atomistic simulations. For this purpose, aluminum, modeled using embedded atom potential of [71], is chosen as the metal of interest and the atomistic simulations are performed using LAMMPS [72]. Ovito [73] is used to visualize the results. The MD box geometry and setup is shown in Fig. 18. A simulation box with dimensions (L x , L y , L z ) = (20.0, 42.2, 1.2) nm is selected. The simulation box has free edges in x-direction and the system is periodic in z-direction. The loading (shown by u x and u y ) is applied on the upper edge in y-direction while the lower one is kept fixed. The grain boundary type (misorientation) is the well-known Σ 5 GB, and the lattice orientation for each of the grains is shown in Fig. 18. The initial crack (L c = 10.67 nm) is introduced in the simulation box which represents an initial defect in the interface to facilitate crack growth at the GB. Further studies may be performed for more realistic GBs, considering the influence of different defects (see [74]). The chosen system orientation inhibits dislocation activity, thus neglecting the need for a plasticity model for regions inside the grain. Coupling the current

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836 15

Fig. 16. Influence of the parameter L x while keeping other interface parameters constant (see Table 1).

Fig. 17. Influence of the parameter L x for Hl = 0.0 (no healing) while keeping other interface parameters constant (see Table 1).

Fig. 18. Left: Geometry and boundary conditions of the MD simulation box. The brown regions are used to highlight the regions, where the normal gap is measured. Analogously, the yellow regions are used for the shear gap. Right: Dimensions of the MD body and the Σ 5 grain boundary. The color coding is based on common neighbor analysis (CNA) implemented in Ovito. Green represents a perfect fcc structure, red represents hcp (i.e. stacking faults) and blue represents the unknown atoms (i.e. dislocation cores and disordered grain boundaries or free surfaces such as crack surface). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

16 S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836

Fig. 19. Front view of the MD simulation under mode I loading. The color coding is the same as in Fig. 18. Table 2 Parameters for the interface geometry. Interface initial crack length L c Interface length L x

Unit

Value

[nm] [nm]

10.0 20.0, 32.0

method with a plasticity model for the grains, such as dislocation density-based formulations [75,76] or phase field models of dislocations [77], represents work in progress. Tractions are measured close to the interface to eliminate the influence of atomic irregularity at the crack surface in calculating the stresses (see brown regions in Fig. 18). These regions are chosen close enough to the GB so that the elastic deformation within the brown region is negligible. The system’s energy is minimized using the over damped method (Sheppard et al. 2008). The internal stress and the grain boundary structure is dynamically relaxed during 10,000 MD steps (step size of 2 fs) under NPT ensemble. The relaxation is performed under zero external stress in the periodic and no load in other directions while the temperature is kept constant at 300 K. The system is then quenched to 20 K under another 10,000 NPT steps with the same conditions as above. Damping factors of 0.1 ps and 0.2 ps are chosen for the thermostat and barostat. After the initialization and relaxation of the simulation box, the system is loaded by displacing the atoms at the ˚ top edge and fixing them at the bottom. The load is applied with a velocity of 0.05 A/ps and the simulation roughly lasts 100,000 steps. During the loading, the velocity rescale method is applied to keep the temperature constant at 20 K. Some snapshots of the system are shown in Figs. 19 and 20, respectively. Intergranular fracture is observed for mode I loading (see Fig. 19). The normal gap here is measured according to gn = gn f − gn0 . On the other hand, for mode II loading (Fig. 20), the two grains slide against each other until they reach the full separation. Interestingly, as the grains slide, new bondings are formed between the atoms which are now getting closer to each other. The shear gap gs (as shown in Fig. 20) is measured based on the average difference of the position of a small selected region of atoms (shown in yellow) inside the brown regions in Fig. 15. Choosing the brown region for measuring the shear gap yield the same results as the yellow region. The same geometry and boundary conditions as the one in the MD simulation are used for the FE calculation (see Fig. 21 and Table 2). Note that to include the influence of the initial crack, some of the interface elements are removed (see the middle part of Fig. 21). As it is shown later on, during the mode II loading, it could be that the leftmost upper nodes become close to the rightmost lower nodes and form new bondings. Therefore, the additional interface elements are considered (shown in blue color in the right-hand side of Fig. 21). The parameters for the interface model are fitted to the results of the MD simulations for mode I and mode II loadings (see Table 3). The results of this comparison are shown in Fig. 22. For mode I loading, the initial slope, the maximum traction, and the final normal gap opening are captured accurately. By only considering the results of mode I loading, the parameters ∆0 , T0 and ∆ f are determined. The

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836 17

Fig. 20. Front view of the MD simulation under mode II. The color coding is the same as in Fig. 18. Table 3 Fitted parameters for the interface model based on atomistic calculations. Maximum traction density T0 Initial interface opening ∆0 Critical interface opening ∆ f Healing length parameter Hl

Unit

Value

[GPa/nm2 ] [nm] [nm] [nm]

2.0 1.2 2.4 1.2

remaining parameter is the healing distance Hl which can be obtained by comparing the results of mode II loading. At this point, more accurate optimization tools may be used to obtain better fitting with more accurate numbers. According to Fig. 20, as the upper grain slides against the lower one, eventually the initial crack heals by forming new bonding. The latter observation explains the constant shear traction in Fig. 22(b). Moreover, as it is shown in Fig. 20, the interface between the GB does not get separated perfectly at the end. This explains the difference between the results of the MD simulation and the FE one at the end of mode II loading (compare blue and black curve in Fig. 22(b)). It should be emphasized that the results reported in Fig. 22 clearly show the anisotropic nature of the GB in normal and shear direction, although a very simple formulation with isotropic damage variable is used for traction density (see 4 parameters reported in Table 1). This suggests that what is measured as an anisotropic behavior is a combination of several structural features at lower scales. Using the current method, not only the correct anisotropic behavior is captured accurately, but also the balance of angular momentum is satisfied. Please note that these MD simulations which are also reported for the first time take several hours (with parallel computing) to be completed

18 S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836

Fig. 21. Finite element setup for comparison with MD simulation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 22. Comparison of the results from molecular dynamic (MD) and finite element (FE) simulation for mode I (φ = 90◦ ) and mode II (φ = 0◦ ) loading. The interface parameters are according to Table 3.

on RWTH-University clusters. The new method accurately predicts the GB behavior in a few seconds on a single CPU pc. To make sure that grain boundary sliding is represented accurately, the system under mode II is unloaded as shown in Fig. 23. Interestingly, the predictions of the model match perfectly with the MD results regarding the remaining shear gap and also the behavior of the GB through the full unloading. Here, no additional internal variable for plasticity nor any additional yield criterion or evolution law is used (see also studies by [78]). For further comparison, the results of the mixed-mode simulations from MD are compared to the FE predictions using the parameters reported in Table 3. The results for the predicted tractions are in good qualitative agreement with MD calculations (see Figs. 24 and 25). For similar studies (mixed-mode loading) using MD simulations, see also [59,79]. Finally, to address the influence of the interface length, the interface size L x is increased while keeping the initial crack length L c unchanged. Using the same parameters for the interface, the results from MD and FE calculations are compared (see Fig. 26). Comparing the MD simulation results with the FE prediction, it is possible to capture the interface behavior using the same data provided in Table 3. A very important point is the behavior of the shear traction in mode II loading which is completely size-dependent. Not only the peak for the shear traction increases by increasing L x , but also the final opening (where shear traction vanishes) is also increasing. The latter observation again can be explained based on the healing effects. The longer the interface length is, the higher is the final gap at which the shear traction vanishes.

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836 19

Fig. 23. Comparison of the results from molecular dynamic (MD) and finite element (FE) simulation for mode II (φ = 0◦ ) loading and unloading. The interface parameters are according to Table 3.

Fig. 24. Comparison of the normal traction from MD and FE simulation for mixed mode loading.

Fig. 25. Comparison of the shear traction from MD and FE simulation for mixed mode loading.

20 S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836

Fig. 26. Comparison of the results from molecular dynamic (MD) and finite element (FE) simulation for different L x . For the interface parameter see Table 3.

Fig. 27. Boundary conditions and meshing of a sample FE model.

6.4. Volume element simulation A simple finite element simulation using the introduced interface model is performed and its results are discussed in this section. A 2D grain morphology and its discretization are shown in Fig. 27. The computations are based on the plane strain assumption. For simplicity, hexagonal grain shapes are considered and discretized via quadrilateral solid elements. Such and similar morphologies show up in many typical applications (such as study of fracture in bone [80]). The interface between the grains (grain boundaries) is enriched by the new interface element formulation. It should be noted that using the current nonlocal method, the bandwidth of the global stiffness matrix increases. In general, for reducing the computational cost it is better to connect the nodes which are and will remain in the vicinity of each other. For simplicity of mesh generation, the nodes in Fig. 27 are connected to all the segments (which is unnecessary). The FE box is loaded in x-direction and other boundary conditions are shown at the left hand side of Fig. 27. A pre-crack is introduced at the top edge to localize the stresses in the middle of the specimen.

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836 21

Fig. 28. Results of the finite element simulation using calibrated interface model.

For the chosen example, the grain size is small enough (about 20 nm) such that different deformation mechanisms are activated at the grain boundaries. Furthermore, in order to focus on the influence of the GB, a simple isotropic elastic model is used for the grain interior (i.e. solid elements). In this simulation, properties of Aluminum are considered for the elastic grains i.e. E = 70 GPa and ν = 0.2. In general, there are also other complexities which have to be taken into account in future works. Among those, one should mention the competition between intergranular and transgranular fracture [81] and interactions of the GB with the neighboring grains through dislocation movements [82–85]. The results of the above simulation in terms of stress contours are presented in Fig. 28. Different rows in Fig. 28 belong to the evolution of different stress components during the loading. The crack in this example tends to propagate in a zig-zag path until it reaches the lower edge. As expected, the maximum stress values are near the crack tip. At the triple junctions (where three grains meet each other at one corner), stress concentration is predicted employing the current model which is also comparable with similar works [86–89]. More investigation on the latter point should be considered in future work. To study different active mechanisms at the GB, similar grain morphology is used but with different boundary conditions. According to Fig. 29, the right edge is under prescribed displacement in x and y directions. The displacement component u y is almost twice the displacement component u x , deforming the system in shear mode mainly. Interestingly, one can observe that although an initial crack is introduced at the top edge, another vertical crack from the lower edge starts to propagate. Moreover, based on the current boundary conditions, cracks tend to develop in (roughly) 45◦ angle. Thanks to the current formulation, one is able to capture GB decohesion together with severe GB sliding as shown in Fig. 29. Further quantitative studies on the obtained reaction force versus the applied displacement are considered as future investigations. It is worth to mention that these computations using

22 S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836

Fig. 29. Results of the finite element simulation using calibrated interface model.

the new nonlocal interface method can still compete with the local methods regarding the computational time. The reason is the very simple formulation of the traction density vector which is very easy to compute. 7. Discussion According to our investigation on GB using MD simulations, nonlocality (for certain applications such as GB) is not a matter of choice but has to be considered to capture the physics at the interface correctly. One important issue is to avoid violation of any basic balance law. Our goal is, first of all, to report about this fact and secondly to introduce a simple approach to capture that in our numerical model. Nevertheless, there are still some open questions and also many possible future developments which are addressed in the current section. 7.1. Comparison with other interface models As mentioned before, models to describe the interface mechanics have been developed to a high level in the literature through the years. Here, we would like to just mention some of the major contributions we found and compare them to our current formulation. • According to Figs. 11–14, the behavior of the obtained traction–separation relation using the current formulation can be tuned easily based on the introduced input parameters. As a result, qualitatively similar behavior can be obtained compared to early cohesive zone models, i.e. the model by Xu and Needleman [10], Tvergaard and Hutchinson [11] as well as Geubelle and Baylor [12]. As an example, by reducing the healing distance Hl , the behavior of the shear traction changes from the bilinear shape [12] to the exponential form [10]. Interestingly, according to Tvergaard and Hutchinson [90], plasticity enhances the interface toughness especially when it comes to mode II loading. These investigations are also in agreement with our results for mode II loading. Compare the results reported in Fig. 22(b) with classical trapezoidal traction–separation relation [90,91]. Please note that these models generally do not satisfy rotational equilibrium as shown by Vossen et al. [24]. Moreover, the interface behavior during unloading or reloading should be addressed properly.

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836 23

• The PPR model developed by Park et al. [14] is a general potential-based interface model which characterizes different fracture energies in each fracture mode. Further numerical comparisons with PPR model would be also interesting especially regarding the work-of-separation in different loading cases (aligned with [25]). Nevertheless, it is worth to mention that in the PPR approach four physical fracture parameters are employed in each fracture mode [92]. On the other hand, we showed that the number of input parameters can be significantly decreased once more information on the interface at smaller scales is taken into account. In that sense, the proposed methodology is also comparable to multi-scaling approach for modeling interface mechanics [93,94]. • One of the cohesive zone models at the presence of large deformation is developed by Van den Bosch et al. [22]. They managed to include mode-dependent work-of-separation into their model. Since in their model the traction vector is co-linear with the gap vector, balance of angular momentum is also satisfied. Nevertheless, the latter assumption is very restrictive for a general interface. Within the current work, we show that this problem can be solved by considering the nonlocal nature of the traction vector at the interface (see also Fig. 7). • A nonlocal cohesive zone was also developed by Paggi and Wriggers [95]. In contrast to the current work, they considered nonlocal interaction perpendicular to the interface. Considering the atomic interaction at smaller scales, this extension can be important in some applications. However, in their work, an anisotropic relation for traction–separation relation is assumed which in general violates the balance of angular momentum if further not further assumptions are considered. It is worth to mention that, the violation of the balance of angular momentum is not due to the poor discretization approach rather weak or simple assumptions at the modeling level. As a consequence, even the node-to-segment cohesive zone developed by the same authors may not satisfy the balance of angular momentum [29]. • A recent approach to tackle the problems regarding interface modeling is to write the free energy of a general interface as a function of the displacement jump (the gap vector) and other internal variables [5,17,86,96]. Later on, to assure thermodynamic consistency of the model, proper evolution laws should be introduced for different mechanisms such as plasticity and damage. The interesting aspect about the current work is to represent plasticity by a combination of damage and healing processes at the grain boundary. Doing so, in contrast to previously developed models, there is no need for evolution laws regarding the plastic displacement jumps. As a consequence, fundamental principles such as the balance of angular momentum and material frame indifference are satisfied (see [17]). At this point, numerical comparisons with models developed by Ottosen et al. [17,20] would be very interesting especially when it comes to prediction of the plastic gap upon unloading the interface. • An important and perhaps essential extension for the description of a general interface at finite displacement jumps is to consider additional tractions related to membrane-like forces (out-of-plane shear forces) acting within the interface (see e.g. models by Ottosen et al. [26], Javili et al. [97] and Heitbreder et al. [27]). It has been shown that considering the aforementioned extension, one is able to solve the issues with the classical cohesive zone modeling (e.g. violation of the balance of angular momentum and material frame indifference). Perhaps, there are also similarities between obtaining the traction via traction density (current nonlocal approach) and considering these additional out-of-plane shear forces. It might be that this effect is already embedded in our formulation. But this remains to be shown. So, it would be interesting to carry out a numerical comparison between this and the presented approach. With respect to the reviewed models, the following new contributions are made. A novel nonlocal interface model is proposed which is rooted in the physical picture of the system. As a result, the model is able to resolve a number of issues of the classical phenomenological cohesive zone models (similar to models developed by Ottosen et al. [26] and Heitbreder et al. [27]). When it comes to grain boundary modeling, on the other hand, the size-dependent results reported in Fig. 26 are of great importance which can be captured by the new proposed model. Therefore, studies on failure initiation and evolution in polycrystalline materials can be improved [98–100]. 7.2. Further future developments In the previous section, we compared the new nonlocal method with previously developed models in the literature. Nevertheless, numerical comparisons between different interface models should be also investigated in future works. Further future developments are as follow:

24 S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836

• In order to model arbitrary crack propagation in the entire body, the proposed methodology has to be extended. The latter point can be seen as an interesting extension and improvement for the nonlocal models such as gradient-extended continuum damage models or phase-field models for fracture [35,101–103]. As a result, one is able to accurately simulate an arbitrary crack propagation within polycrystalline materials [104,105]. Moreover, the interaction of intragranular plasticity with the GB has to be taken into account [82,106–108]. • The principles of the current work are also similar to other nonlocal techniques developed in the literature such as peridynamics [33,109,110]. Further comparisons and investigations in this direction are also insightful. • Other discretization methods can be utilized for the proposed method. Even using a meshless method should be possible. In this way, different nodes from opposite faces of the interface are connected depending on the nature of the problem and influence area. • As discussed in Section 4.2, the idea of traction density or nonlocal interactions at the interface can be considered in different applications even at the macroscale (it is not just limited to small scales). As an example, one can mention delamination in composites and study of fiber bridging [61,62]. The definition of the traction density has to adopted based on the application. This adoption depends on the physics of the interface at hand (see also [63]). • The proper thermo-mechanical coupling will be another interesting and also necessary aspect [111]. This can also be achieved by performing different MD simulations to study temperature influences on the interface mechanical behavior. Furthermore, one may also include the influence of friction within the interface model [100]. • Study of fatigue and cyclic loading under more complex boundary conditions and different type of morphologies is also an important point which results in optimum design of materials for specific application [112]. 8. Conclusion A new method is proposed to model the mechanics of an interface. By introducing a quantity called traction density, it is possible to capture the complex behavior of an arbitrary interface. The traction–separation relation is obtained by integrating the traction density along the interface. The formulation of the traction density is relatively simple and requires physical knowledge about the interface under study. By applying the proposed formulation it is possible to reduce the number of required parameters (inputs) to model the interface. Furthermore, setting up complicated phenomenological models can be avoided. The proposed nonlocal formulation is easy to implement and conventional cohesive zone models can also be extended to include these features, if it is necessary. Finally, satisfying basic balance laws at finite displacement jumps can be assured. The latter point is a well-known problem for anisotropic interface models. In this work, the focus is on the grain boundaries. Previous and current studies using molecular dynamics simulations show two main active phenomena at the grain boundary: intergranular fracture and grain boundary sliding. The coupling of these two phenomena is a challenging task which has been addressed in the current work. Current MD simulations also suggest that when it comes to grain boundary sliding, size-dependent results should be considered. In other words, the obtained shear traction for the sliding depends on the interface length. In the current formulation, it is possible to capture this for the very first time using the healing agent. Interestingly enough, using only two MD simulations in mode I and II loading, the interface formulation can be calibrated. Predictions obtained by using the calibrated interface model for the mixed-mode simulation are in a very good agreement with the MD simulations. The latter point shows how well the new model can capture the behavior of the grain boundary in a consistent way. Acknowledgments Financial support of Subproject A6 of the Transregional Collaborative Research Center SFB/TRR 87 by the German Research Foundation (DFG) is gratefully acknowledged. Appendix A. Convergence study The results for convergence study with respect to number of elements at the interface n el and number of integration points for each element n I P are shown in Figs. 30 and 31, respectively.

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836 25

Fig. 30. Convergence study with respect to number of elements n el keeping the number of integration points constant n I P = 13. The interface parameters are according to Table 1.

Fig. 31. Convergence study with respect to number of integration points n I P keeping the number of interface elements constant n el = 20. The interface parameters are according to Table 1.

By increasing the number of elements at the interface, the obtained tractions converge to a final solution. The healing effect is the main reason behind the oscillatory behavior of the shear traction for very rough mesh size. Similar behavior and explanations are also valid for the convergence study with respect to the number of integration points (see Fig. 31). Appendix B. Balance of linear and angular momentum Considering the assumptions in Eqs. (10) and (11), the balance of linear ( p) and angular momentum (L) resulted from the traction forces at the crack interface is satisfied as ∫ ∫ p = ∫Γ + ∫t + dΓ + + Γ − t − dΓ ∫− ∫ = ∫Γ + ∫Γ − T − dΓ − dΓ + + Γ − Γ + T + dΓ + dΓ − = ∫Γ + ∫Γ − (T − + T + ) dΓ − dΓ + (27) = Γ + Γ − (T − T ) dΓ − dΓ + = 0,

26 S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836

and L

= = = = =

∫ ∫ + + + − − − ∫Γ + ∫x × +t dΓ− + −Γ − x + × t∫ dΓ ∫ + − + − ∫Γ + ∫Γ − x + × T − dΓ dΓ − + Γ+− Γ + x × T dΓ dΓ T dΓ dΓ ∫Γ + ∫Γ − (x − x ) × − G × T dΓ dΓ + Γ+ Γ− 0,

(28)

since G and T are collinear. References [1] S. Osovski, A. Needleman, A. Srivastava, Intergranular fracture prediction and microstructure design, Int. J. Fract. (2019). [2] A. Akbari, P. Kerfriden, S. Bordas, On the effect of grains interface parameters on the macroscopic properties of polycrystalline materials, Comput. Struct. 196 (2018) 355–368. [3] D. Farkas, S. Van Petegem, P. Derlet, H. Van Swygenhoven, Dislocation activity and nano-void formation near crack tips in nanocrystalline Ni, Acta Mater. 53 (2005) 3115–3123. [4] J.S.K.-L. Gibson, S. Rezaei, H. Rueß, M. Hans, D. Music, S. Wulfinghoff, J.M. Schneider, S. Reese, S. Korte-Kerzel, From quantum to continuum mechanics: studying the fracture toughness of transition metal nitrides and oxynitrides, Mater. Res. Lett. 6 (2) (2018) 142–151. [5] S. Rezaei, D. Jaworek, J.R. Mianroodi, S. Wulfinghoff, S. Reese, Atomistically motivated interface model to account for coupled plasticity and damage at grain boundaries, J. Mech. Phys. Solids 124 (2019) 325–349. [6] P.P. Camanho, C.G. Davila, M.F. de Moura, Numerical simulation of mixed-mode progressive delamination in composite materials, J. Compos. Mater. 37 (16) (2003) 1415–1438. [7] S. Rezaei, M. Arghavani, S. Wulfinghoff, N.C. Kruppe, T. Brögelmann, S. Reese, K. Bobzin, A novel approach for the prediction of deformation and fracture in hard coatings: Comparison of numerical modeling and nanoindentation tests, Mech. Mater. 117 (2018) 192–201. [8] S. Rezaei, S. Wulfinghoff, S. Reese, Prediction of fracture and damage in micro/nano coating systems using cohesive zone elements, Int. J. Solids Struct. 121 (2017) 62–74. [9] I. Khisamitov, G. Meschke, Variational approach to interface element modeling of brittle fracture propagation, Comput. Methods Appl. Mech. Engrg. 328 (2018) 452–476. [10] X.-P. Xu, A. Needleman, Numerical simulations of fast crack growth in brittle solids, J. Mech. Phys. Solids 42 (9) (1994) 1397–1434. [11] V. Tvergaard, J.W. Hutchinson, On the toughness of ductile adhesive joints, J. Mech. Phys. Solids 44 (5) (1996) 789–800. [12] P. H.G.eubelle, J. S.B.aylor, Impact-induced delamination of composites: a 2D simulation, Composites B 29 (5) (1998) 589–602. [13] M. Ortiz, A. Pandolfi, Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis, Internat. J. Numer. Methods Engrg. 44 (9) (1999) 1267–1282. [14] K. Park, G.H. Paulino, J.R. Roesler, A unified potential-based cohesive model of mixed-mode fracture, J. Mech. Phys. Solids 57 (6) (2009) 891–908. [15] F. Confalonieri, U. Perego, A new framework for the formulation and validation of cohesive mixed-mode delamination models, Int. J. Solids Struct. (2019). [16] P. Carrara, L.D. Lorenzis, A coupled damage-plasticity model for the cyclic behavior of shear-loaded interfaces, J. Mech. Phys. Solids 85 (2015) 33–53. [17] N.S. Ottosen, M. Ristinmaa, Thermodynamically based fictitious crack/interface model for general normal and shear loading, Int. J. Solids Struct. 50 (22) (2013) 3555–3561. [18] R. Dimitri, M. Trullo, L.D. Lorenzis, G. Zavarise, Coupled cohesive zone models for mixed-mode fracture: A comparative study, Eng. Fract. Mech. 148 (2015) 145–179. [19] M. Leuschner, F. Fritzen, J. van Dommelen, J. Hoefnagels, Potential-based constitutive models for cohesive interfaces: Theory, implementation and examples, Composites B 68 (2015) 38–50. [20] N.S. Ottosen, M. Ristinmaa, J. Mosler, Fundamental physical principles and cohesive zone models at finite displacements – Limitations and possibilities, Int. J. Solids Struct. 53 (2015) 70–79. [21] T. Heitbreder, N.S. Ottosen, M. Ristinmaa, J. Mosler, Consistent elastoplastic cohesive zone model at finite deformations – Variational formulation, Int. J. Solids Struct. 106–107 (2017) 284–293. [22] M. van den Bosch, P. Schreurs, M. Geers, Identification and characterization of delamination in polymer coated metal sheet, J. Mech. Phys. Solids 56 (11) (2008) 3259–3276. [23] O. van der Sluis, B. Vossen, J. Neggers, A. Ruybalid, K. Chockalingam, R. Peerlings, J. Hoefnagels, J. Remmers, V. Kouznetsova, P. Schreurs, M. Geers, Advances in delamination modeling of metal/polymer systems: Continuum aspects, in: J.E. Morris (Ed.), Nanopackaging: Nanotechnologies and Electronics Packaging, Springer International Publishing, Cham, 2018, pp. 83–128. [24] B. Vossen, P. Schreurs, O. van der Sluis, M. Geers, On the lack of rotational equilibrium in cohesive zone elements, Comput. Methods Appl. Mech. Engrg. 254 (2013) 146–153. [25] J. Mosler, I. Scheider, A thermodynamically and variationally consistent class of damage-type cohesive models, J. Mech. Phys. Solids 59 (8) (2011) 1647–1668. [26] N.S. Ottosen, M. Ristinmaa, J. Mosler, Framework for non-coherent interface models at finite displacement jumps and finite strains, J. Mech. Phys. Solids 90 (2016) 124–141.

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836 27 [27] T. Heitbreder, N.S. Ottosen, M. Ristinmaa, J. Mosler, On damage modeling of material interfaces: Numerical implementation and computational homogenization, Comput. Methods Appl. Mech. Engrg. 337 (2018) 1–27. [28] F. Parrinello, G. Marannano, G. Borino, A thermodynamically consistent cohesive-frictional interface model for mixed mode delamination, Eng. Fract. Mech. 153 (2016) 61–79. [29] M. Paggi, P. Wriggers, Node-to-segment and node-to-surface interface finite elements for fracture mechanics, Comput. Methods Appl. Mech. Engrg. 300 (2016) 540–560. [30] G. Pijaudier-Cabot, Z.P. Bažant, Nonlocal damage theory, J. Eng. Mech. 113 (10) (1987) 1512–1533. [31] Z.P. Bažant, M. Jirásek, Nonlocal integral formulations of plasticity and damage: Survey of progress, J. Eng. Mech. 128 (11) (2002) 1119–1149. [32] R.H.J. Peerlings, R. De Borst, W.A.M. Brekelmans, J.H.P. De Vree, Gradient enhanced damage for quasi-brittle materials, Internat. J. Numer. Methods Engrg. 39 (19) (1996) 3391–3403. [33] S. Silling, E. Askari, A meshfree method based on the peridynamic model of solid mechanics, Comput. Struct. 83 (17) (2005) 1526–1535. [34] A.C. Hansen-Dörr, R. de Borst, P. Hennig, M. Kästner, Phase-field modelling of interface failure in brittle materials, Comput. Methods Appl. Mech. Engrg. 346 (2019) 25–42. [35] L. Sharma, R.H.J. Peerlings, P. Shanthraj, F. Roters, M.G.D. Geers, An FFT-based spectral solver for interface decohesion modelling using a gradient damage approach, Comput. Mech. (2019). [36] J. Nordmann, K. Naumenko, H. Altenbach, A damage mechanics based cohesive zone model with damage gradient extension for creep-fatigue-interaction, in: Advances in Engineering Plasticity and Its Application IX, in: Key Engineering Materials, vol. 794, Trans Tech Publications, 2019, pp. 253–259. [37] V.P. Nguyen, C.T. Nguyen, S. Bordas, A. Heidarpour, Modelling interfacial cracking with non-matching cohesive interface elements, Comput. Mech. 58 (5) (2016) 731–746. [38] R. Dimitri, L. De Lorenzis, P. Wriggers, G. Zavarise, Nurbs- and t-spline-based isogeometric cohesive zone modeling of interface debonding, Comput. Mech. 54 (2) (2014) 369–388. [39] R. Janisch, N. Ahmed, A. Hartmaier, Ab initio tensile tests of al bulk crystals and grain boundaries: Universality of mechanical behavior, Phys. Rev. B 81 (2010) 184108. [40] X. Zhou, N. Moody, R. Jones, J. Zimmerman, E. Reedy, Molecular-dynamics-based cohesive zone law for brittle interfacial fracture under mixed loading conditions: Effects of elastic constant mismatch, Acta Mater. 57 (16) (2009) 4671–4686. [41] S. Lee, V. Sundararaghavan, Calibration of nanocrystal grain boundary model based on polycrystal plasticity using molecular dynamics simulations, Int. J. Multiscale Comput. Eng. 8 (5) (2010) 509–522. [42] L. Guin, J.L. Raphanel, J.W. Kysar, Atomistically derived cohesive zone model of intergranular fracture in polycrystalline graphene, J. Appl. Phys. 119 (2016) 245107. [43] S.P. Patil, Y. Heider, C.A.H. Padilla, E.R. Cruz-Chú, B. Markert, A comparative molecular dynamics-phase-field modeling approach to brittle fracture, Comput. Methods Appl. Mech. Engrg. 312 (2016) 117–129. [44] M.G. Elkhateeb, Y.C. Shin, Molecular dynamics-based cohesive zone representation of ti6al4v/tic composite interface, Mater. Des. 155 (2018) 161–169. [45] J.J. Möller, E. Bitzek, R. Janisch, H. ul Hassan, A. Hartmaier, Fracture ab initio: A force-based scaling law for atomistically informed continuum models, J. Mater. Res. 33 (22) (2018) 3750–3761. [46] D. Warner, F. Sansoz, J. Molinari, Atomistic based continuum investigation of plastic deformation in nanocrystalline copper, Int. J. Plast. 22 (4) (2006) 754–774. [47] H. Song, Y. Li, Atomic simulations of deformation mechanisms of crystalline Mg/amorphous Mg–Al nanocomposites, Phys. Lett. A 379 (36) (2015) 2087–2091. [48] A. Elzas, T. Klaver, B. Thijsse, Cohesive laws for shearing of iron/precipitate interfaces, Comput. Mater. Sci. 152 (2018) 417–429. [49] J. Gong, A.J. Wilkinson, Sample size effects on grain boundary sliding, Scr. Mater. 114 (2016) 17–20. [50] H. Conrad, L. Rice, The cohesion of previously fractured fcc metals in ultrahigh vacuum, Metall. Trans. 1 (11) (1970) 3019–3029. [51] A.A. Shirzadi, H. Assadi, E.R. Wallach, Interface evolution and bond strength when diffusion bonding materials with stable oxide films, Surf. Interface Anal. 31 (2001) 609–618. [52] K. Khaledi, S. Rezaei, S. Wulfinghoff, S. Reese, A microscale finite element model for joining of metals by large plastic deformations, C. R. Mecanique 346 (8) (2018) 743–755, computational modeling of material forming processes Simulation numérique des procédés de mise en forme. [53] K. Khaledi, S. Rezaei, S. Wulfinghoff, S. Reese, Modeling of joining by plastic deformation using a bonding interface finite element, Int. J. Solids Struct. 160 (2018) 68–79. [54] J. Mergheim, P. Steinmann, Phenomenological modelling of self-healing polymers based on integrated healing agents, Comput. Mech. 52 (3) (2013) 681–692. [55] A.A. Alsheghri, R.K.A. Al-Rub, Thermodynamic-based cohesive zone healing model for self-healing materials, Mech. Res. Commun. 70 (2015) 102–113. [56] S.A. Ponnusami, J. Krishnasamy, S. Turteltaub, S. van der Zwaag, A cohesive-zone crack healing model for self-healing materials, Int. J. Solids Struct. 134 (2018) 249–263. [57] S. Forest, Micromorphic approach for gradient elasticity, viscoplasticity, and damage, J. Eng. Mech. 135 (3) (2009) 117–131. [58] T. Brepols, S. Wulfinghoff, S. Reese, Gradient-extended two-surface damage-plasticity: Micromorphic formulation and numerical aspects, Int. J. Plast. 97 (2017) 64–106. [59] A. Elzas, B. Thijsse, Cohesive laws describing the interface behaviour of iron/precipitate interfaces under mixed loading conditions, Mech. Mater. 129 (2019) 265–278.

28 S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836 [60] Y. Pan, F. Tian, Z. Zhong, A continuum damage-healing model of healing agents based self-healing materials, Int. J. Damage Mech. 27 (5) (2018) 754–778. [61] Y. Li, S. Reese, J.-W. Simon, Modeling the fiber bridging effect in cracked wood and paperboard using a cohesive zone model, Eng. Fract. Mech. 196 (2018) 83–97. [62] D. Höwer, B.A. Lerch, B.A. Bednarcyk, E.J. Pineda, S. Reese, J.-W. Simon, Cohesive zone modeling for mode I facesheet to core delamination of sandwich panels accounting for fiber bridging, Compos. Struct. 183 (2018) 568–581, in honor of Prof. Y. Narita. [63] B.F. Sørensen, E.K. Gamstedt, R.C. Østergaard, S. Goutianos, Micromechanical model of cross-over fibre bridging – prediction of mixed mode bridging laws, Mech. Mater. 40 (4) (2008) 220–234. [64] B.F. Sørensen, T.K. Jacobsen, Large-scale bridging in composites: R-curves and bridging laws, Composites A 29 (11) (1998) 1443–1451. [65] J. Reinoso, M. Paggi, A consistent interface element formulation for geometrical and material nonlinearities, Comput. Mech. 54 (6) (2014) 1569–1581. [66] S.A. Silling, M. Epton, O. Weckner, J. Xu, E. Askari, Peridynamic states and constitutive modeling, J. Elasticity 88 (2) (2007) 151–184. [67] A. Javili, R. Morasata, E. Oterkus, S. Oterkus, Peridynamics review, Math. Mech. Solids 0 (0) (2018) 1081286518803411. [68] N. Zhu, D.D. Meo, E. Oterkus, Modelling of granular fracture in polycrystalline materials using ordinary state-based peridynamics, Materials 9 (2016) 997. [69] S. Gur, M.R. Sadat, G.N. Frantziskonis, S. Bringuier, L. Zhang, K. Muralidharan, The effect of grain-size on fracture of polycrystalline silicon carbide: A multiscale analysis using a molecular dynamics-peridynamics framework, Comput. Mater. Sci. 159 (2019) 341–348. [70] S. Jafarzadeh, Z. Chen, F. Bobaru, Peridynamic modeling of intergranular corrosion damage, J. Electrochem. Soc. 165 (2018) C362–C374. [71] Y. Mishin, D. Farkas, M.J. Mehl, D.A. Papaconstantopoulos, Interatomic potentials for monoatomic metals from experimental data and ab initio calculations, Phys. Rev. B 59 (1999) 3393–3407. [72] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117 (1) (1995) 1–19. [73] A. Stukowski, Visualization and analysis of atomistic simulation data with ovito–the open visualization tool, Modelling Simulation Mater. Sci. Eng. 18 (1) (2010) 015012. [74] A. Tehranchi, W. Curtin, Atomistic study of hydrogen embrittlement of grain boundaries in nickel: I. Fracture, J. Mech. Phys. Solids 101 (2017) 150–165. [75] F. Roters, P. Eisenlohr, L. Hantcherli, D.D. Tjahjanto, T.R. Bieler, D. Raabe, Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: Theory, experiments, applications, Acta Mater. 58 (4) (2010) 1152–1211. [76] J. Rezaei Mianroodi, R. Peerlings, B. Svendsen, Strongly non-local modelling of dislocation transport and pile-up, Phil. Mag. (2016) 1–17. [77] J.R. Mianroodi, B. Svendsen, Atomistically determined phase-field modeling of dislocation dissociation, stacking fault formation, dislocation slip, and reactions in fcc systems, J. Mech. Phys. Solids 77 (2015) 109–122. [78] Z.H. Aitken, D. Jang, C.R. Weinberger, J.R. Greer, Grain boundary sliding in aluminum nano-bi-crystals deformed at room temperature, Small 10 (1) (2014) 100–108. [79] B. Paliwal, M. Cherkaoui, An improved atomistic simulation based mixed-mode cohesive zone law considering non-planar crack growth, Int. J. Solids Struct. 50 (2013) 3346–3360. [80] L. Lin, X. Wang, X. Zeng, An improved interfacial bonding model for material interface modeling, Eng. Fract. Mech. 169 (2017) 276–291. [81] D.M. Bond, M.A. Zikry, Differentiating between intergranular and transgranular fracture in polycrystalline aggregates, J. Mater. Sci. 53 (8) (2018) 5786–5798. [82] P. van Beers, G. McShane, V. Kouznetsova, M. Geers, Grain boundary interface mechanics in strain gradient crystal plasticity, J. Mech. Phys. Solids 61 (12) (2013) 2659–2679. [83] S. Wulfinghoff, A generalized cohesive zone model and a grain boundary yield criterion for gradient plasticity derived from surfaceand interface-related arguments, Int. J. Plast. 92 (2017) 57–78. [84] F. Bormann, R. Peerlings, M. Geers, B. Svendsen, A computational approach towards modelling dislocation transmission across phase boundaries, Phil. Mag. 99 (2019) 2126–2151. [85] E. Alabort, D. Barba, S. Sulzer, M. Lißner, N. Petrinic, R. Reed, Grain boundary properties of a nickel-based superalloy: Characterisation and modelling, Acta Mater. 151 (2018) 377–394. [86] C.F. Dahlberg, J. Faleskog, C.F. Niordson, B.N. Legarth, A deformation mechanism map for polycrystals modeled using strain gradient plasticity and interfaces that slide and separate, Int. J. Plast. 43 (2013) 177–195. [87] T.-L. Cheng, Y.-H. Wen, J.A. Hawk, Diffuse interface approach to modeling crystal plasticity with accommodation of grain boundary sliding, Int. J. Plast. 114 (2019) 106–125. [88] C. Pu, Y. Gao, Y. Wang, T.-L. Sham, Diffusion-coupled cohesive interface simulations of stress corrosion intergranular cracking in polycrystalline materials, Acta Mater. 136 (2017) 21–31. [89] A. Simone, C.A. Duarte, E. Van der Giessen, A generalized finite element method for polycrystals with discontinuous grain boundaries, Internat. J. Numer. Methods Engrg. 67 (8) (2006) 1122–1145. [90] V. Tvergaard, J.W. Hutchinson, The influence of plasticity on mixed mode interface toughness, J. Mech. Phys. Solids 41 (6) (1993) 1119–1135. [91] I. Scheider, W. Brocks, Simulation of cup–cone fracture using the cohesive model, Eng. Fract. Mech. 70 (14) (2003) 1943–1961.

S. Rezaei, J. Rezaei Mianroodi, K. Khaledi et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112836 29 [92] K. Park, G.H. Paulino, Cohesive zone models: A critical review of traction-separation relationships across fracture surfaces, Appl. Mech. Rev. 64 (2011). [93] B. Vossen, P. Schreurs, O. van der Sluis, M. Geers, Multi-scale modeling of delamination through fibrillation, J. Mech. Phys. Solids 66 (2014) 117–132. [94] K. Matouš, M.G. Kulkarni, P.H. Geubelle, Multiscale cohesive failure modeling of heterogeneous adhesives, J. Mech. Phys. Solids 56 (4) (2008) 1511–1533. [95] M. Paggi, P. Wriggers, A nonlocal cohesive zone model for finite thickness interfaces – Part I: Mathematical formulation and validation with molecular dynamics, Comput. Mater. Sci. 50 (5) (2011) 1625–1633. [96] M.E. Gurtin, L. Anand, Nanocrystalline grain boundaries that slip and separate: A gradient theory that accounts for grain-boundary stress and conditions at a triple-junction, J. Mech. Phys. Solids 56 (1) (2008) 184–199. [97] A. Javili, P. Steinmann, J. Mosler, Micro-to-macro transition accounting for general imperfect interfaces, Comput. Methods Appl. Mech. Engrg. 317 (2017) 274–317. [98] H.D. Espinosa, P.D. Zavattieri, A grain level model for the study of failure initiation and evolution in polycrystalline brittle materials. Part I: Theory and numerical implementation, Mech. Mater. 35 (3) (2003) 333–364. [99] T. Luther, C. Könke, Polycrystal models for the analysis of intergranular crack growth in metallic materials, Eng. Fract. Mech. 76 (15) (2009) 2332–2343. [100] I. Benedetti, M. Aliabadi, A three-dimensional cohesive-frictional grain-boundary micromechanical model for intergranular degradation and failure in polycrystalline materials, Comput. Methods Appl. Mech. Engrg. 265 (2013) 36–62. [101] F. Cazes, A. Simatos, M. Coret, A. Combescure, A cohesive zone model which is energetically equivalent to a gradient-enhanced coupled damage-plasticity model, Eur. J. Mech. A Solids 29 (6) (2010) 976–989. [102] C.V. Verhoosel, R. de Borst, A phase-field model for cohesive fracture, Internat. J. Numer. Methods Engrg. 96 (1) (2013) 43–62. [103] A.C. Hansen-Dörr, L. Wilkens, A. Croy, A. Dianat, G. Cuniberti, M. Kästner, Combined molecular dynamics and phase-field modelling of crack propagation in defective graphene, Comput. Mater. Sci. 163 (2019) 117–126. [104] J. Clayton, J. Knap, Phase field modeling of directional fracture in anisotropic polycrystals, Comput. Mater. Sci. 98 (2015) 158–169. [105] T.-T. Nguyen, J. Réthoré, J. Yvonnet, M.-C. Baietto, Multi-phase-field modeling of anisotropic crack propagation for polycrystalline materials, Comput. Mech. 60 (2) (2017) 289–314. [106] P. van Beers, V. Kouznetsova, M. Geers, Defect redistribution within a continuum grain boundary plasticity model, J. Mech. Phys. Solids 83 (2015) 243–262. [107] N.A. Giang, A. Seupel, M. Kuna, G. Hütter, Dislocation pile-up and cleavage: effects of strain gradient plasticity on micro-crack initiation in ferritic steel, Int. J. Fract. 214 (2018) 1–15. [108] F. Bormann, R. Peerlings, M. Geers, On the competition between dislocation transmission and crack nucleation at phase boundaries, Eur. J. Mech. A Solids 78 (2019) 103842. [109] E. Askari, F. Bobaru, R.B. Lehoucq, M.L. Parks, S.A. Silling, O. Weckner, Peridynamics for multiscale materials modeling, J. Phys. Conf. Ser. 125 (2008) 012078. [110] D. De Meo, N. Zhu, E. Oterkus, Peridynamic modeling of granular fracture in polycrystalline materials, Trans. ASME, J. Eng. Mater. Technol. 138 (4) (2016). [111] M. Fagerström, R. Larsson, A thermo-mechanical cohesive zone formulation for ductile fracture, J. Mech. Phys. Solids 56 (10) (2008) 3037–3058. [112] S. Roth, G. Hütter, M. Kuna, Simulation of fatigue crack growth with a cyclic cohesive zone model, Int. J. Fract. 188 (1) (2014) 23–45.