20-Multiscale Methods For A Fracture A Review

J. Indian Inst. Sci.
A Multidisciplinary Reviews Journal

ISSN: 0970-4140 Coden-JIISAD
© Indian Institute of Science 2017.
ARTICLE
REVIEW
Multiscale Methods for Fracture: A Review⋆
P. R. Budarapu1,2* and T. Rabczuk3,4
Abstract | The global response of a system is often governed by the

material behaviour at smaller length scales. Investigating the system
mechanics at the smallest scale does not always provide the complete
picture. Therefore, in the ambitious objective to derive the overall full-
scale global response using a bottom-up approach, multiscale methods
coupling disparate length and time scales have been evolved in the past
two decades. The major objective of the multiscale methods is to reduce
the computational costs by coupling the inexpensive coarse-scale/con-
tinuum based models with expensive fine-scale models. The fine-scale
region is employed in the critical areas, such as crack tips or core of the
dislocation. To improve the efficiency the fine-scale domain is adaptively
adjusted as the defects propagate. As a result, the accuracy of the fine-
scale model is combined with the efficiency of the coarse-scale model,
arriving at a computationally efficient and accurate multiscale model.
Currently, multiscale methods are applied to study problems in numer-
ous fields, involving multiphysics. In this article, we present an overview
of the multiscale methods for fracture applications. We discussed the
techniques to model the coarse- and fine-scale domains, details of the
coupling methods, adaptivity, and efficient coarse-graining techniques.
The article is concluded with comments on recent trends and future
scope.
Keywords: Multiscale methods, Multiphysics, Computational fracture, Atomistic simulations, Coarse
graining, Adaptivity
1 Introduction to reveal the fundamental mechanics of material

The global response of a system is often governed failure at nano scales by modeling the atom-to-
by the material behaviour at smaller length scales. atom interactions, due to their small dimensions 1
School of Mechanical
For example, the macroscopic properties of a of the order of angstroms (Å), they are still pro- Sciences, Indian Institute
material such as toughness, strength, ductility, hibitively expensive to be employed in industrial of Technology Bhubaneswar,
thermal and electrical conductivity, and chemical applications2,3. A plausible alternative to reduce Bhubaneswar 752050,
India. 2 Department
diffusion are strongly influenced by defects like the computational demand is to couple the con- of Aerospace Engineering,
cracks and dislocations, which are initiated and tinuum scale with the discrete scale using a mul- Indian Institute of Science,
evolved at the nano scales. In the ambitious objec- tiscale approach. In such paradigms, defects are Bangalore 560012, India.
3
Institute of Structural
tive to derive the overall full-scale global response explicitly modeled at the sub-scales, whilst a self- Mechanics, Bauhaus
using a bottom-up approach, the sub-scale behav- consistent continuum model elsewhere. Several University of Weimar,
iour has to be accurately computed. Therefore, numerical models dealing with multiple spatial 99423 Weimar, Germany.
4
Division of Computational
understanding the phenomena of material failure and temporal length scales have been proposed Mechanics, Faculty of Civil
across multiple length scales has been the major in the past two decades1,4–12. Most of the coupling Engineering, Ton Duc Thang
research focus in the material science and engi- methods and simulations are focused on models University, Ho Chi Minh City,
⋆
Vietnam. Dedicated to the
neering community for many years1. Although of intact materials (without cracks). The trans- researchers of IISc.
molecular dynamics (MD) simulations promise fer of information through different length scales * pattabhib@gmail.com
J. Indian Inst. Sci. | VOL 97:3 | 339–376 September 2017| journal.iisc.ernet.in 13

P. R. Budarapu and T. Rabczuk
for problems involving material failure and finite phenomenological multiscale model, where the
temperatures still remains a challenging task. In elastic properties are computed using direct
this paper, an overview of the multiscale methods homogenization and subsequently evolved using
for fracture applications is presented. We denote a simple three-parameter orthotropic continuum
the sub-scale domain the “fine scale region” and damage model. A unified regularization scheme
the continuum domain the “coarse scale region”. was employed in the context of constitutive law
Multiscale methods are initiated through the rescaling and the staggered nonlocal approach.
quasi-continuum method (QCM)1,13, by directly A hierarchical higher order multiscale cohesive
coupling the fine-scale region to the coarse-scale zone model (MCZM) is introduced in23, to simu-
region. In the QCM, the continuum degrees of late the fracture in crystalline solid, using up to
freedom need to be exactly located at the posi- third-order Cauchy–Born rules and Barycentric
tions of the atoms at the interface, which is finite element method to construct shape func-
achieved by fine grading of the continuum mesh tions for hexagonal shaped cohesive zones.
around the coupling region. The QCM has also Lawrimore et al.22 studied the mechanical
been very successful at linking two continuum responses of low volume fraction Polyvinyl Alco-
scales, for example, for fibrous materials14 and hol /Montmorillonite nano-composites using
is readily capable of including quantum effects a hierarchical multiscale modeling method,
through density functional theory (QCDFT)15. by bridging the MD with finite element analy-
Beex et al.16 have investigated four variants of the sis (FEA). MD computations of interfaces were
quasi-continuum method for their use in planar used to calibrate a traction–separation relation,
beam lattices which can also experience out-of which was then upscaled based on a cohesive
plane deformation. Different frameworks are zone model (CZM)57–59. Paggi et al.20 estimated
compared to the direct lattice computations for the influence of micro-cracking and power-loss
three truly multiscale test cases in which a single in photovoltaic modules based on a global–local
lattice defect is present in an otherwise perfectly multi-physics multi-scale approach. The micro-
regular beam lattice. The virtual-power-based scale damage mechanisms, particularly matrix/
quasi-continuum method is adopted for lattice inter-phase fracture and fiber sliding in brittle
models in which bond failure and subsequent ceramics, are studied in60 based on multiscale
frictional fiber sliding are incorporated, which methods. Nguyen et al.50 have presented a com-
are of significant importance for fibrous materials putational homogenization procedure for cohe-
such as paper, cardboard, textile, and electronic sive and adhesive crack modeling of materials
textile17. Bond failure and fiber sliding are non- with a heterogeneous micro structure, for crack
local dissipative mechanisms, which are treated propagation under cyclic loading with numerical
with a mixed formulation in which the kinematic analysis of the convergence characteristics of the
variables as well as the internal history variables multiscale method and treatment of macroscopic
are interpolated. Multiscale methods can be cat- snapback in a multiscale simulation. Greco et al.61
egorized into hierarchical18–23, semi-concur- have proposed a concurrent multiscale method
rent24–31 and concurrent methods4,6–12,32–43, see to overcome the existing limitations on homog-
Fig.1. enization in the presence of strain localization in
In hierarchical multiscale methods, see Fig. 1a, masonry structures. They adopted a multilevel
the information is passed from the fine-scale domain decomposition approach equipped with
to the coarse-scale; but not vice versa. Com- an adaptive zooming-in criterion for detecting
putational homogenization31,44–52 is a classical the zones affected by strain localizations.
up-scaling technique. Hierarchical multiscale However, hierarchical multiscale methods
approaches are very efficient. Therefore, hierar- based on computational homogenization are not
chical methods have been successfully applied to well suited to model fracture. One basic assump-
study various problems, ranging from multiphase tion of homogenization theories is the existence of
flow in porous media53–55 to polymer nano-com- disparate length scales62: LCr ≪ LRVE ≪ LSpec ,
posites (PNC). An iterative multiscale finite vol- where LCr, LRVE and LSpec are the size of the:
ume method for the simulation of multiphase crack, representative volume element (RVE)
flow in fractured porous media in the context of a and specimen, respectively. The first condition
hierarchical fracture modeling framework is dis- is violated for problems involving fracture, as
cussed in18. A hierarchical multiscale approach is LCr is of the order of LRVE. Furthermore, peri-
employed in21 to characterize the material behav- odic boundary conditions (PBC), often specified
iour of the heat-affected zone (HAZ) in welded at the fine-scale, cannot be used when a crack
connections. Liu et al.56 developed a regularized touches a boundary as the displacement jump
340
13 J. Indian Inst. Sci.| VOL 97:3 | 339–376 September 2017 | journal.iisc.ernet.in
Multiscale Methods for Fracture: A Review
(a) Nano (c)

(b)
F
Micro scale
(RVE)
C Continuum
Micro C P
Crack front
Macro F B
A
C F P
P 4 3
5 2
F 1 A
Macro scale 6 7 Interface Handshake
Nano scale (VAC) (BSM) (BDM)
coarse fine coarse fine coarse fine
Figure 1: Schematics of a hierarchical, b semi-concurrent, and c concurrent multiscale methods. In the

hierarchical methods the information exchange happens only from fine scale to coarse scale, whereas
the interaction is two way in case of semi-concurrent and concurrent multiscale methods. Note a definite
region of coupling in the concurrent multiscale methods, which does not exist in the semi-concurrent mul-
tiscale methods.
in that boundary violates the periodic boundary element software package ABAQUS. A macro-
conditions. micro model has been implemented in ANSYS
Figure 1b illustrates the basic idea of semi- software67, to predict matrix cracking evolution
concurrent multiscale methods, where the in laminates under in-plane loading by treating
information is passed from the fine-scale to the transverse cracks as separate discontinuities
the coarse-scale and vice versa. A classical semi- in the micro-model. Kerfriden et al.68 discussed
concurrent multiscale method is the FE 224–27 a technique to reduce the computational burden
method, originally developed for intact materials associated with the simulation of localized failure
and later extended to problems involving mate- in global–local framework. Recently, Ojo et al.69
rial failure9,45,48,49,63. Oliver et al.64 discussed the have proposed a non-local adaptive discrete
computational strategies to affordably solve mul- empirical interpolation method combined with
tiscale fracture problems using FE 2 approach. modified hp-refinement for order reduction of
Computational efficiency of the semi-concurrent molecular dynamics systems.
multiscale methods is similar to concurrent mul- Numerous concurrent multiscale methods
tiscale methods65. The key advantage of semi- have been developed that can be classified into
concurrent multiscale methods over concurrent ‘Interface’ coupling methods and ‘Handshake’
multiscale methods is their flexibility, i.e., their coupling methods, see Fig. 1c. The coupling hap-
ability to couple two different software packages, pens along an interface in case of interface cou-
for example, MD software to FE software66. Zhu pling methods32, whereas, a definite region of
et al.30 introduced a nonlinear semi-concurrent coupling exists for the handshake coupling meth-
multiscale method to model crack propagation ods7. Interface coupling methods are not efficient
evolving from micro-structure for non-linear for dynamic applications, as avoiding spuri-
material behaviour based on an asymptotic ous wave reflections at the ’artificial’ interface is
expansion homogenization combined with the problematic. On the other hand, since the coarse
semi-concurrent finite element approach. Modi- region does not exist everywhere, adaptive adjust-
fied periodic boundary conditions and spheri- ment of the fine-scale region as the defects propa-
cal grains’ generation procedure are devised for gate in the handshake method is cumbersome.
nonlinear material model with post-failure stage. Some of the concurrent multiscale methods have
A semiconcurrent multiscale computational been extended to modeling fracture12,39,70. Talebi
homogenization method for the simulation of et al.43 have developed a concurrent coupling
hydro-mechanical problem for quasi-brittle scheme coupling molecular dynamics to extended
materials is proposed in31, coupling MD-mod- finite element method through the bridging
els with continuum models. Silani et al.29 have domain method, to model three-dimensional
developed a semi-concurrent multiscale method cracks and dislocations at the atomistic level.
to estimate the pre-localized damage initia- Budarpu et al.70 have proposed a solid shell-based
tion and propagation in the fully exfoliated clay/ concurrently coupled three-dimensional adap-
epoxy nanocomposite, where the methodology tive multiscale method (3DAMM) to simulate
has been implemented in the commercial finite complex crack growth patterns in thin-walled
341
J. Indian Inst. Sci. | VOL 97:3 | 339–376 September 2017 | journal.iisc.ernet.in 13
structures. The material in the bulk is modeled possesses four nearest neighbors. Therefore,
using a hybrid solid shell formulation relying two Silicon atoms share four outer most atoms
on the combined use of the enhanced assumed to form covalent bonds. When the Sun light
strain (EAS) and the assumed natural strain equal to or more than a photon is incident on a
(ANS) methods71–75. The authors developed a PV cell, electrons from the outermost orbit are
computational framework in MATLAB by trig- excited and released from their regular orbits,
gering Large-scale Atomic/Molecular Massively creating holes and free electrons. The free elec-
Parallel Simulator (LAMMPS)76 through sys- trons are free to move throughout the crystal77.
tem command. A photovoltage is generated by creating a poten-
In this paper, we present a review of the tial barrier across the moving charges. On the
recent trends in multiscale methods for fracture other hand, micro-cracks in Silicon cells can lead
applications. The article is structured as follows: to power losses up to 21%78–82. Therefore, the
multiscale methods are introduced in Sect. 1. Var- mechanics of a PV cell involves mechanical, ther-
ious techniques to model fine- and coarse-scale mal and electrical fields, requiring multiphysics-
domains are discussed in Sect. 2 and 3, respec- based techniques for accurate analysis20,83–85.
tively. Section 4 is dedicated to the techniques on Furthermore, MD simulations of multiphys-
coupling the coarse- and fine-scale domains of ics models demand potential functions such as
concurrently coupled multiscale methods. Crack charge optimized many body (COMB) poten-
nucleation/growth criteria and the techniques tial86, considering the combined effect of various
for adaptive adjustment of the fine-scale region fields.
are discussed Sect. 5. Computer implementation The total potential energy (E tot(q, r)), consid-
steps of a three-dimensional enhanced bridging ering the effects of charge transfer, is expressed
scale method based concurrently coupled mul- in terms of the electrostatic energies (E es(q, r)),
tiscale method, in the MATLAB frame work are short-range interactions (E short(q, r)), van der
presented in Sect. 6. The article is concluded in Waals interactions (E vdW(r)), and correction
Sect. 7 with a discussion on perspective future terms (E corr(r)), where q and r represent the
developments of multiscale methods. charges and atom positions, respectively86:
E tot (q, r) = E es (q, r) + E short (q, r)

2 Fine‑scale modeling (1)
In this section, various techniques to closely + E vdW (r) + E corr (r).
analyze the mechanics of fracture, particularly
The application of digital filters to split the
around the crack tip at lower scales are sum-
energy spectrum of an atomistic zone simulated
marized. Popular techniques include atomistic
with molecular dynamics into low- and high-
models, virtual atom cluster models and repre-
energy components is discussed in87. The effect of
sentative volume element approaches.
the grain orientation on the fracture behavior of
polycrystalline silicon in micro-electro-mechan-
2.1 Atomistic models ical systems has been investigated in88, based on
The basic structure of solid materials at nano a multiscale model combining the discontinu-
scales can be obtained by the periodic arrange- ous Galerkin method and extrinsic cohesive law
ment of the unit lattice. Atoms in the lattice struc- describing the fracture process.
ture are bonded together by the van der Waals Techniques to simulate the fracture using
forces of attraction. Atomistic models are parti- atomistic models under static and dynamic sce-
cle-based techniques, where the mechanics are narios are summarized in this section. Static con-
simulated based on the atom–atom interactions. ditions are achieved by minimizing the system
Furthermore, material characteristics depend on potential energy, whereas velocity Verlet scheme89
the arrangement of atoms in the crystal lattice is a popular technique to estimate the atom veloc-
and the forces of interaction such as mechanical, ities and positions.
electrical, chemical, and thermal forces. There-
fore, various types of potential functions are
2.1.1 Molecular dynamics
required to model the atom–atom interactions of
In molecular dynamics, the objective is to deter-
various materials.
mine the atom positions rα (t), velocities vα (t),
Consider a material such as Silicon used in
and their accelerations aα (t), for the given initial
the photovoltaic (PV) applications. Atoms in
conditions. Each atom is assumed to be a classical
the Silicon unit cell are arranged in the diamond
particle obeying Newton’s laws of mechanics. The
cubic lattice structure, where each Silicon atom
342
governing equations are derived from the Lagran- system potential energy will be minimum when
gian equations mentioned below: the first derivative of the potential function with
respect to the positions of the atoms goes to zero.
d ∂L ∂L
− = 0, α = 1, 2, 3, . . . , nA , (2) Therefore, the internal and external forces act-
dt ∂ ṙα ∂rα ing on an atom after the minimization of the
where nA represents the total number of atoms. potential energy are given by12:
The Lagrangian of the fine-scale domain can be
nA nA
estimated as follows: 1 ∂V (rαβ ) ∂rαβ
Fint = −
A A
2 ∂rαβ ∂r
n n α=1 β� =
mα ṙ2 α
(9a)
L= + Uα (r), (3) n A
2 ∂V (rαβ ) rα − rβ
α=1 α=1 =− and
∂rαβ rαβ
where mα is the mass of the atoms and Uα (r) is β� =α
the potential energy of the fine scale, which can
∂W ext
be estimated as sum of all the bond potentials: Fext = − , (9b)
  ∂r
nA nA nA
� � 1 � respectively. The residual forces on each atom are
Uα = φα = V (rαβ ), (4) estimated as R = Fint − Fext. Details of deriva-

2
α=1 α=1 β� =α
tion of Eq. (9) and the computer implementation
where φα is the potential energy associated with steps are explained in.12
atom α and β’s denote all the neighbors of atom Molecular dynamics simulations became
α. Therefore, substituting Eq. (3) in Eq. (2) yields important tools to understand the key mecha-
the equations of motion in Newtonian form: nisms in wide range of applications. There-
fore, over the years MD simulations attracted
∂U(rα )
mα r̈α = = Fα , α = 1, 2, 3, . . . , nA , (5) researchers in several fields, starting from early
∂rα simulations of Abraham et. al.2,90 for engineer-
where Fα is the internal force vector acting on ing applications to modern applications such as
atom α. Considering the initial conditions, Eq. (5) estimating the material properties91–99, simu-
is solved for the trajectories of the atomic motion late problems involving multiphysics20,100–105,
in the current configuration. Brief Verlet algo- model the mechanics of polymer nano compos-
rithm steps include the following: (1) solve Eq. ites22,106–109 and bio-medical applications110–114,
(6) for the current atom positions to name a few. Therefore, several exclusive soft-
wares have been evolved in the recent years to
1
rα (t + �t) = rα (t) + vα (t)�t + aα (t)�t 2 , (6) carry out MD simulations in various fields. Some
2 of the popular softwares include (1) the open
where (2) the atom velocites are estimated from source Large-scale Atomic/Molecular Massively
Parallel Simulator (LAMMPS)76, with included
1
(7)
vα (t + �t) = vα (t) + [aα (t) + aα (t + �t)]�t, potential functions to simulate the mechancis of
2 soft and solid-state materials and coarse-grained
in which (3) the accelerations in the current time systems. (2) Accelerated molecular dynamics
step aα (t + �t) in Eq. (7) are calculated from the simulations (ACEMD)115,116: an open source
interaction potential function and Eq. (5). There- software to perform molecular dynamics simula-
fore, knowing the accelerations, atom velocities tions of proteins, oligosaccharides, nucleic acids,
and positions can be estimated using Eqs. (7) and and synthetic polymers consisting of Chem-
(6), respectively. istry at HARvard Macromolecular Mechanics
(CHARMM), Amber forcefields and can run
on NVIDIA graphics processing units (GPUs)
2.1.2 Molecular statics
optimized with CUDA. CUDA is a parallel com-
In molecular statics (MS), the aim is to determine
puting platform and programming model that
the positions of the atoms for the given boundary
makes using a GPU for general purpose comput-
conditions, by minimizing the system potential
ing simple and elegant. (3) Abalone117, a general
energy, expressed as
purpose open source molecular modeling pro-
= W int − W ext , (8) gram focused on molecular dynamics of bio-
polymers, which can also interact with external
where W int
represents the internal energy of the
quantum chemical programs ORCA, NWChem,
system and W ext is the external contribution. The
CP2K, and PC GAMESS/Firefly. (4) GROningen
343
MOlecular Simulation (GROMOS)110 is a soft- function φρ for a triangular lattice (see Fig. 2) can
ware for dynamic modeling of bio-molecules. be defined as12:
On the other hand, GROningen MAchine for
7 7
Chemical Simulation (GROMACS)118 is a par- φVAC 1 V (r1β ) 1
φρ = = √ = φ1β (10)
allel message-passing implementation of a MD V0 2 2
3a /2 2
β=2 β=2
program which is useful for bio(macro)molecules
in aqueous environment. (5) AMBER119 is a open Therefore, using the definition of φρ from Eq.
source package of computer programs for apply- (10), the internal nodal forces can be expressed as
ing molecular mechanics, normal mode analysis, follows:12
molecular dynamics, and free energy calculations
to simulate the structural and energetic proper- ∂φρG ∂u
ties of molecules. Apart from the above simu- FIint ≈ − wG
∂u ∂uIC
G
lation softwares, there exist several supporting
7
MD packages, for example, (1) PACKMOL120 is ∂φρG ∂uC α
a package for building initial configurations for ≈− wG (11)
∂uαC ∂uIC
G α=1
molecular dynamics simulations and (2) Visual
G
molecular dynamics (VMD)121 is a popular visu- n
7
∂φρG
alization software to post process the MD simula- =− wG NI (Xα ).
∂uαC
tion results. G=1 α=1
∂φ
The term ∂uCρin Eq. (11) can be evaluated for
α
each atom α in the VAC as given below:
2.2 Virtual atom cluster model
α=1
Virtual atom cluster (VAC) model is based on
considering the symmetry of a crystal structure, ∂φρ ∂φ12 r12i ∂φ13 r13i ∂φ14 r14i
where a cluster of atoms is taken as the represent- C
= + +
∂u1i ∂r r ∂r r ∂r14 r14
ative model of the whole lattice structure122,123. 12 12 13 13
(12)
As a result, all the calculations can be performed ∂φ15 r15i ∂φ16 r16i ∂φ17 r17i
+ + +
with reference to the representative cluster instead ∂r15 r15 ∂r16 r16 ∂r17 r17
of the whole lattice, which leads to improved
α = 2–7
computational efficiency. Since the locations of
atoms in the cluster do not represent the exact ∂φρ ∂φ1α r1αi
=− , (13)
locations of the atoms, the representative cluster ∂uC
αi
∂r1α r1α
is called a virtual atom cluster (VAC). The same
where i is the index of the coordinate axes. The
inter atomic potential as in the full-scale atom- ∂φ
detailed derivation of the term ∂uCρ is given in
istic model can be used in the VAC model as 12 α
appendix of . Knowing the internal nodal forces
well12,122,124. A full-scale atomistic model will be
in Eq. (11), the minimization problem can be
realized when the VAC assumes the structure of
solved for the coarse-scale solution by minimiz-
the underlying lattice. Therefore, VAC is an effi-
ing the potential energy for the given bound-
cient coarse-graining technique to improve the
ary conditions. An extension of the VAC-based
computational efficiency.
coarse-graining scheme to study the dynamic
A schematic of VAC-based coarse-scale model
fracture through a multiscale model is developed
in two dimensions is shown in Fig. 2. The total
in.124
potential energy of a fine-scale system as shown
in Fig. 2a is given by the sum of all bond poten-
tials φα, estimated using Eq. (4). Consider an 2.3 Representative volume element
equivalent coarse-scale model based on the VAC, approach
illustrated in Fig. 2b. Since the fine-scale and Representative volume element approach is pop-
coarse-scale models are equivalent, their potential ularly used in the semi-concurrent multiscale
energy must be equal. This is achieved by defin- methods (see Fig. 1b) to bridge the meso-scale
ing a distributed energy density function φρ122. to the macro-scale29. For fracture analysis based
Considering the periodic nature of the lat- on the micro-mechanical approach using RVE,
tice, φρ is defined as the potential energy of a VAC a damage parameter is estimated in meso-scale,
divided by the volume of the VAC. For a homoge- which will be sent back to the macro-scale for
neous lattice, each VAC consists of a single atom further analysis. The boundary conditions for the
and its volume is that of the unit cell of the lat- RVE are extracted from the macro-scale model.
tice. Therefore, the distributed energy density Huang et al.125 have analyzed the strain hardening
344
and multiple-cracking tensile fracture behavior where E(j) denotes the ensemble average based
of engineered cementitious composites, based on on j realizations, E(2j) indicates the ensemble
a multiscale method. The authors used a multi- average obtained using twice the number of reali-
linear crack bridging relationship at a lower scale zations, and tol2 is a saturation tolerance. The
based on analytical crack bridging analysis for a convergence can be guaranteed for 4 realizations
single crack, whereas a representative volume ele- with a convergence error of less than 1%.
ment model was used at the upper scale using the
extended finite element method. A global–local
multiscale finite element method is employed 3 Coarse‑scale model
in126 to study the interaction of nanotubes and The coarse-scale domain might be discretized
matrix at the nanoscale, based on building a 3D with classical techniques like the finite ele-
finite element model of a representative volume ment method, the Partition of Unity Finite
element around the crack tip. Paggi et al.20 devel- Element Method (PUFEM)127,128, meshfree meth-
oped a multi-physics and multi-scale approach to ods6,129–132 or partition-of-unity enriched meth-
study micro-cracking and power-loss in photo- ods such as the extended Finite Element Method
voltaic modules. (XFEM)133–141, the Smooth Finite Element
The size of the RVE depends on the dimension Method (SFEM)138,142–144, the Generalized Finite
of the heterogeneities (d), which are expected to Element Method (GFEM)145–150, the extended
be much smaller than the dimension of the RVE. Element Free Galerkin method (XEFG)151–158,
Therefore, the suitable size of the RVE for a par- the Cracking Particles Method159–165, the Phan-
ticular problem can be arrived after few succes- tom node method (PNM)166–172 or the Numerical
sive entanglement tests. For example, consider Manifold Method (NMM)173,174, the phase-field
the ensemble average of the elastic modulus esti- methods102,175–178, Peridynamics179,180, to name
mated based on different spatially random sam- a few, apart from the isogeometric analysis with
ples using different increasing RVE sizes, until the high-order approximation techniques51,142,181–186
below stagnation condition is satisfied29: with exact geometry.
Alternative techniques to model fracture
(1)
|�E�l ′ − �E�l |
(1) in shells include a FEM-based computational
(1)
< tol1, (14) method for the fracture of plates and shells on
�E�l the basis of edge rotation and load control as
where E(1) described in187. Rabczuk et al.163 have developed a
l is the ensemble average for an RVE
of size l, �E�(1) meshfree method for thin shells with finite strains
l ′ is the ensemble average for the RVE
size of l ′, and tol1 is a stagnation tolerance. A and arbitrary evolving cracks, eliminating the
reasonable number of realizations for the estima- membrane locking using a cubic or fourth-order
tion of the ensemble average can be estimated polynomial basis. However, third-order com-
based on the below saturation criterion29: plete formulations lead to the large support sizes,
increasing the computational cost188. An extrin-
|�E�(2j) − �E�(j) | sic basis to increase the order of completeness
< tol2, (15)
�E�(2j) of the approximation reduces the support size at
4 3
5 2
1
6 7
VAC
VAC Atom
Neighbour atom
4 3
5 2 Coarse scale element
6
1
7 Gauss point
(a) (b) (c)
Figure 2: A demonstration of VAC coarse-scale model in two dimensions. a Atomistic model with trian-
gular lattice as on the (111) plane of an fcc material. b Equivalent continuum model with the VAC being
placed at a particular Gauss point. c Details of the VAC.
345
the cost of adding more degrees of freedom per based on the partition of unity concept127,145 and
node. Nonetheless, numerical experiments show through additional nodal parameters for the ele-
a reduced CPU time for several problems. On the ments cut by the crack. The central idea of XFEM
other hand, linear dependence of the shape func- is to decompose the displacement field into a
tions deteriorates the conditioning of the final continuous part uC and a discontinuous part uD:
system of equations188. Reinoso et al.74 devised
and implemented a 7-parameter shell element u(X) = uC (X) + uD (X) (16)
for geometrically nonlinear analysis of layered The discontinuous part (uD)
is estimated intro-
CFRP composites. Nguyen et al.189 proposed ducing the additional information into the FE
an extended isogeometric element formulation interpolation through the local partition of unity
(XIGA) based on the Kirchhoff–Love theory, approach by adding an enrichment, whereas the
for the analysis of through-the-thickness cracks continuous part (uC) is the standard FE interpo-
in thin shell structures based on Non-Uniform lation. The approximation of the displacement
Rational B-Splines (NURBS). field for nC cracks with nT crack tips reads
Verhoosel et al.63 employed a partition of C
n
unity-based cohesive zone finite element model
u(X) = NI (X)uI + NI (X)ψIK (X)aIK
to mimic crack nucleation and propagation in a
I∈S K =1 I∈S C
piezoelectric continuum, through a multiscale
framework to appropriately represent the influ- nT
NP

M
ence of the microstructure on the response of a + NI (X) φPI (X)bM
PI ,
miniaturized component. Plews et al.190 have M=1 I∈S T P=1 (17)
proposed a generalized finite element approach where S is the set of nodes in the entire discre-
for predicting localized, nonlinear, thermoplastic tization, S C is the set of nodes associated with
behavior and residual stresses and deformations completely cracked elements, S T is the set of
in structures subjected to intense heating. A mul- nodes around the crack tip, NI is the standard
tiscale reduction technique to describe shale gas shape functions, ψIK (X) is the enrichment func-
transport in fractured media is discussed in191. T is the enrichment function
tion of Kth crack, φPI
The matrix is described by upscaled models and for the crack tip P, aI , and bPI are the additional
the interaction between the matrix and the frac- degrees of freedom to be solved. Further details
tures is modeled through the generalized multi- on XFEM and its applications can be found in
scale finite element method192. the excellent review papers by Belytschko et al.195,
In this section, popular techniques to model Fries and Belytschko196, Karihaloo and Xiao197,
cracks in the continuum, such as extended finite Rabczuk66 or the book by Mohammadi198.
element method, meshless methods, phase field
method and the phantom node method are dis-
3.2 Extended Meshless methods
cussed. Mathematical formulation of the last
Meshless methods (MM) evolved after the intro-
three methods are presented in detail.
duction of Element Free Galarkin (EFG) method
by Belytschko130. They are based on the idea of
3.1 Extended finite element method using the nodes/particles in the zone of influence
Finite element methods, a powerful class of tech- of a selected point to construct its approximation
niques to study a wide range of problems, fail to space. The particles are not connected; hence, it is
model the problems involving strong and weak particularly easy to simulate complex phenomena
discontinuities with a jump in the displacement like fracture129,157,161,165,199–202. Popular meshless
and strain field, because of their smooth interpo- methods include Meshless Local Petrov-Galerkin
lation character. In other words, computational (MLPG)203,204, the reproducing kernel particle
failure mechanics involving fracture related to the method (RKPM)131,132, radial point interpolation
initiation and propagation of crack, which falls method (RPIM)205,206, to name a few. A mesh-
under the strong discontinuities, and interface less variational multiscale method for thermo-
problems involving interactions of solid–solid mechanical material failure is discussed in101. The
and fluid–fluid (weak discontinuities) and solid– displacement and temperature fields are enriched
fluid (strong discontinuities) are cumbersome with step-functions and appropriate crack tip
to simulate using the FE techniques. Therefore, enrichment accounting for fine-scale features.
to handle linear and nonlinear crack openings, Yang et al.207 have developed a meshless col-
a very flexible extended finite element method location method based on the differential repro-
has been developed by Belytschko et al.193,194 ducing kernel (DRK) approximation, where the
346
derivative of the shape function of reproduc- interpolant function NIcont (X) is defined ask
ing kernel (RK) approximants is replaced by follows:
a set of differential reproducing conditions to
avoid the complex direct differentiation. On the NIcont (X) = N̂I (X) + N̄I (X), (20)
other hand, in the EFG130,208 method, the shape where N̂ (s) is estimated based on the quartic
functions are directly differentiated. Therefore, spline function. N̄I (x) in Eq. (20), is introduced
the number of degrees of freedom are reduced to impose the nth order reproducing conditions:
in DRK approximation, improving the com-
putational efficiency of the DRKP method as N̄I (X) = wa (X − XI )PT (X − XI )z̄(X), (21)
compared to the RKPM. The Kronecker delta where the weight function wa (X − XI ) is cen-
property is not satisfied by the shape functions of tered at XI , defined by the normalized Gaussian
the RK approximants. To resolve the difficulty in function, and PT (X − XI ) are the nth order poly-
the DRK approximation, a meshless collocation nomial functions. For an nth order complete pol-
method based on the DRKP method has been ynomial basis, a set of reproducing conditions can
introduced by ChingPing et al.104,105,209,210, satis- be obtained to determine z̄i (X), i = 1, 2, . . . , nB ,
fying the Kronecker delta property. where nB are the total number of basis functions
Consider a body in two-dimensional vec- given by (n + 1)(n + 2)/2. Further details of esti-
tor space R2 with boundary Ŵ, as shown in Fig. 3. mation of N̂I (X) and N̄I (X) functions are dis-
The crack is denoted by c on the surface Ŵc. cussed in42.
The enriched shape functions in Eq. (19) are
expressed as the product of the standard shape
3.3 Displacement field
function and the sign function H:
A displacement field that is discontinuous at the
crack(s) and continuous elsewhere in the domain NIenr (X) = NIcont (X)H(fI (X)), (22)
is a proper choice to model fracture in meshless
where
methods. Therefore, the total displacement field
is decomposed into a standard/continuous and

1∀ξ > 0
discontinuous/enriched part: H(ξ ) = (23)
−1∀ξ < 0,
u = ucont + uenr , (18) and fI (X) = n0 · (X − XI ), n0 is the normal of
the split nodes in the initial configuration, esti-
where ucontis the continuous component and
mated from the fine scale.
uenr is the discontinuous component. The mesh-
free approximation of Eq. (18) is given by
3.3.1 Variational formulation
uh (X) = NIcont (X)uI + NIenr (X)qI , The governing equations in weak form can be
c
(19)
I∈S I∈S
stated as follows: find u ∈ U , ∀δu ∈ U0, such that,
where NIcont and NIenr are the displacement inter-
δW = δWint − δWext = 0, (24)
polation/approximation functions in the con-
tinuous and discontinuous domains, respectively, where u = ū on Ŵu and u ∈ H(�), δWint is the
see Fig. 4. uI and qI indicate the nodal parameters first variation of the internal energy, and δWext is
associated with the continuous and discontinu- the virtual work from the external forces. Let the
ous displacement fields, respectively. The DRKP test functions δuh (X) be defined as
(a)
u,
t,
t
(b)
c,
c n+ -
- c
+ +
n- c
Figure 3: Domain and surface descriptions of a body a displacements, tractions and the crack along with
their surfaces b a closeup of the crack tip indicating the directions on the normals on the crack surfaces.
347
3.4 Phase‑field method for fracture

δuh (X) = NIcont (X)δuI + NIenr (X)δqI .
I∈S I∈S c Phase field (PF) approaches have been extensively
(25) used over the years for several different engineer-
Therefore, based on Eqs. (19), (24) and (25), the ing applications175,176,211–213. The PF approach has
discrete system of equations can be obtained as42 been proven to be very effective and accurate to
simulate complex fracture patterns177,178,214–216,
KIJ · dJ = fIext , (26) even at nano-scales217. On the other hand, PF
where KIJ is the stiffness matrix, and indicatefIext model requires very fine discretizations to accu-
the external force vector, and dJ is the vector con- rately capture failure in solids, which is compu-
taining the nodal parameters. The stiffness matrix tationally expensive. Based on the value of the
is estimated as phase field parameter (d), the orientation of the
crack surface and the approximate location of the
uq
KIJuu crack tip can be identified. Giovanardi et al.218

KIJ
KIJ = qu qq d�, (27)
�0 \Ŵ0 c KIJ KIJ proposed a global–local approach by coupling
XFEM in the caorse scale with the phase field
where model at the fine scale. Yingjun et al.217 developed
a multiscale method based on PF approach, to
∂NIcont (X) ∂NJcont (X) simulate the microstructure evolution of materi-
KIJuu = C , (28a)
∂X ∂X als at nanoscales. They simulated the morphol-
ogy of microcrack propagation in single crystal
uq ∂NIcont (X) ∂NJcont (X)H(fJ (X)) materials under tensile strain with a fixed grip
KIJ = C , (28b)
∂X ∂X condition, by coupling phase field crystal with an
external field method.
qu ∂NIcont (X)H(fI (X)) ∂NJcont (X) The central idea of the phase field approach of
KIJ = C , (28c)
∂X ∂X brittle fracture consists of regularizing the sharp
crack topology within a diffusive crack zone of
qq ∂NIcont (X)H(fI (X)) ∂NJcont (X)H(fJ (X)) width l, see Fig. 5, where the regularization of the
KIJ = C ,
∂X ∂X sharp crack representation is depicted. Therefore,
(28d) it is possible to define a scalar-valued function
where C is the matrix of material constants. The that accounts for the stiffness degradation such
expressions for the nodal forces are given by159: that d, with d : loc0 × [0, t] → [0, 1], refer
175,176
for further details.

fIext = NIcont (X)T bd� Within the framework of brittle frac-
�0 \Ŵ0c ture211,214,215, the potential energy function for
(29) a cracked body can be defined as the sum of the
+ NJenr (X)T bd�. deformation strain energy �(E) integrated over
�0 \Ŵ0c ph
the domain 0 and the critical fracture energy Gc
1.2 1
0.8
1
0.6
0.8
0.4
0.6 0.2
0
0.4
−0.2
0.2 −0.4
−0.6
0
−0.8
−0.2 −1
−1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1
Figure 4: Shape functions in a the continuous domain and b the discontinuous domain. Picture repro-
duced with permission from42.
348
Figure 5: a Sharp and diffusive crack modeling. Left: discrete crack discontinuity in the continuum
domain. Right: smeared discontinuity in the continuum domain based on the PF approach. b Diffusive
crack modeling solution for the one-dimensional crack problem.
integrated along the crack path Ŵc, which can be in which

expressed as
G d (u, d, δ d) = −2(1 − d)δ d�(E) d�loc
0
�loc
�(u, d) = g(d)�(E) d�loc0 0

1

�loc
0
+ Gc l dδ d + ∇ X d · ∇ X (δ d) d�loc
0
d2 l2

Gc l �loc

2
+ + |∇ X d| d�loc
0
0
(32)
�loc
0
2 l2 = 0, ∀δ d ∈ Vd ,
(30)
+ �ext (u), where Vu = δu ∈ [H 1 (�loc0 )] : δu = 0 on
where ext is the contribution due to the pre- ∂ B0,u denotes the space of admis-
scribed external actions (from the global level), sible displacement variations, and
and g(d) = [1 − d]2 + K identifies the mono- Vd = {δ d ∈ H1 (B0 ) | δ d = 0 on Ŵc } stands for
tonically decreasing degradation function, K ≈ 0 the space of admissible test functions for the
being a residual stiffness parameter. phase field. The second Piola-Kirchhoff stress
tensor is defined as S := ∂E �.
In the global scale, the interpolation of the
3.4.1 Finite element formulation displacement field (u) along with its variation
The first variation of Eq. (30) with respect to the (δu ) and increment (u) using the standard tri-
independent fields E and u is given by175,176 linear shape functions can be expressed as
G u (u, δu, d) = Gint
u u
− Gext u ≈ Nd, δu ≈ Nδd, and �u ≈ N�d. (33)

∂� ∂E The phase field interpolation (d), its variation
= g(d) : δu d�loc
0 + δ�ext (u)
loc
�0 ∂E ∂u (δ d), and increment (d) at the element level are
(31)
= 0, ∀δu ∈ Vu ,
349
approximated using the same shape functions the continuum domain into two sub-domains
corresponding to the kinematic field: �0 = �0(+) �0(−). Correspondingly, two phan-
p p p
tom domains are defined: �0 = �0(+) �0(−).
d = Nd, δ d = Nδ d, and �d = N�d. (34) Since the elements in the two sub domains do not
The gradient of the phase field (∇X d), its varia- share any nodes in common, their displacements
tion (∇X δ d), and increment (∇X d) are inter- are independent, resulting in the expected discon-
polated through a suitable operator Bd, as tinuity across the cross surface.
mentioned below: Through the definition of f as the signed
distance measured from the crack surface,
∇X d = Bd d, ∇X δ d = Bd δ d, ∇X �d = Bd �d. W0+ , Wp− , W0−, and Wp+ as the nodes belonging
(35) p p
to �0(+) , �0(−) , �0(−), and �0(+), respectively, the
Applying the displacement and the phase field discontinuous interpolation of the displacement
discretization at the local level to Eqs. (31)–(32) field is given by
and performing the consistent linearization of the
residual vectors, a coupled system of equations u(X, t) = uI (t)NI (X)H(f (X))
can be obtained: I∈{W0+ ,Wp− }
u
+ uJ (t)NJ (X)H(−f (X)),

G (u, δu, d) !
= 0. (36)
G d (u, d, δ d) J ∈{W0− ,Wp+ } (37)
where H is the Heaviside function. In line
3.5 Phantom node method
with167,219, the standard approximation of the dis-
In the phantom node method (PNM)166–169,
placements on each part of the cracked element
when an element is completely cut by a crack,
�0(+) and �0(−), which are extended to their cor-
the displacement field can be represented as con- p p
responding phantom domains �0(−) and �0(+)
tinuous on each part of the cracked element and
introduces the continuous displacement field.
discontinuous across the crack surface. There-
The displacement jump between the two flanks
fore, the crack kinematics can be obtained by
of the crack can be computed by taking the dif-
overlapping elements168–172 using the extra nodes
ference of the displacement fields of the two
known as phantom nodes. Therefore, (1) the dis-
domains of the cracked element.
placement field is discontinuous across the crack
but independently continuous on each part of
the cracked element. Hence, the discontinuous 4 Coupling techniques
element is replaced by two elements with the In this section, we present the techniques to cou-
additional phantom nodes, which requires only ple the coarse- and fine-scale domains. Two pop-
a small modification in existing finite element ular techniques, bridging scale method (BSM)
codes; (2) the associated shape functions in a and bridging domain method (BDM) to cou-
cracked element are the same as the shape func- ple the continuum with atomistic domains are
tions of an intact element, and (3) the elements discussed.
adjacent to the cracked elements do not require
any modification. Because of the above advan-
4.1 Bridging scale method
tages, the computer implementation of the phan-
An overview of the multiscale method based on
tom node method is particularly easy.
the bridging scale concept is presented in220 with
Consider an arbitrary continuous body with
an emphasis on complex material systems. In this
a surface of discontinuity Ŵc. According to the
section, an outline of a three-dimensional multi-
phantom node method170, the kinematics of
scale method based on enhanced bridging scale
a cracked element can be described by super-
method to model fracture is presented.
imposing two separate displacement fields,
Consider a three-dimensional multiscale
which are active only in a determined region
model shown in Fig. 7 for the adaptive simula-
of the domain. Consequently, a completely
tion of crack growth. The MD model in Fig. 7a
cut element can be represented as an union:
assumes Silicon in the fine-scale domain. The
elem
0 = elem1
0 ∪ elem2
0 , of two elements sepa-
material of the coarse-scale region is modeled
rated along the crack surface, see Fig. 6a–b. The
based on a solid-shell as shown in Fig. 7b. In the
superscript ‘elem’ refers to the considered phan-
diamond cubic lattice structure of Silicon shown
tom element and ‘elem1’ and ‘elem2’ denotes the
in Fig. 7c, each atom possesses four nearest neigh-
sub elements after splitting. This formalism is
bors. Fig. 7b shows the modeling details of the
expressed by setting that the crack surface divides
coarse-scale region highlighting the estimation
350
of stiffness matrix of a solid shell element. The hence, uαC is sufficient to model the deformation
phantom node method170 can be used to model in the coarse-scale region. On the other hand, in
the crack surfaces in the coarse-scale region. the fine-scale region, both coarse- and fine-scale
Crack originates from the coarse-scale components are required. Let the coarse-scale
region and the crack tip is captured in the fine- (see Sect. 3) displacement uαC of an atom α be rep-
scale region. In Fig. 7a, c, the fine-scale region resented by a set of FEM basis functions defined
is formed by the atomistic model, made up of over a set of nC nodal points,
diamond cubic lattice structure of Silicon. Vari-
nC
ous techniques to model the fine-scale domain
are summarized in Sect. 2. The initial crack in uαC = NI (Xα )uIC , (39)
the fine-scale region is created by deleting the I=1
bonds between the atoms on the crack surface where NI (Xα ) is the shape functions defined at
and updating the neighbor list accordingly. Ghost node I , estimated at the αth atom with the mate-
atoms located on the boundary of the coarse rial coordinate Xα, and uIC is the continuum dis-
region, but within the cutoff radius of the atoms placement vector at node I .
in the fine region (see Fig. 7a), are used to enforce In the bridging scale method, the coupling
the boundary conditions for the fine-scale solu- conditions are realized by enforcing the displace-
tion. A finite-discrete element method combining ment boundary conditions on the ghost atoms,
the advantages of both the finite element method see Fig. 8. The positions of the ghost atoms are
and the discrete element method, coupling by interpolated from the coarse-scale solution. Let β
means of ghost particles, is discussed in221. be the index of the ghost atoms; the correspond-
In the two-scale model, the total displacement ing ghost atom displacements are estimated as
field uα of an atom α is decomposed into coarse-
nC
and fine-scale components:
uβC = NI (Xβ )uIC , (40)
uα = uαC + uαA , (38) I=1
where uαC is the coarse-scale component and uαA is where NI (Xβ ) are the shape functions defined at
the fine-scale component. The fine-scale compo- node I, estimated at the βth atom with material
nent uαA is the difference between the actual posi- coordinates Xβ.
tion of an atom α and the interpolated position The bridging scale method has been applied
of the coarse scale. Therefore, uαA is insignificant to study various physical phenomenon. A meso-
in the regions far away from the crack tip, and scopic bridging scale (MBS) method, multiscale
(a) straight crack p 0

8 7 8* 7* 8 7
4 3 4 3
f(X)>0 4* 3* 2
f(X)=0
5 5 1 6 5*
f(X)<0
6 6*
1 2 1 2 1* 2*
0 0 p
(b) angled crack

p
8 7 8* 7 8 0
7*
4 3 4* 3 4 2 3*
) >0
f( X 1
0
)=
f( X 5 6 5 6 5* 6*
0
X )<
f(
1 2 1 2 1* p
2*
0 0
phantom node real node
Figure 6: Schematic representation of a cracked element using the phantom node method, for a a
straight and b angled cracks.
351
(b) Solid shell element

(a) 8 7
4 3
C
Coarse scale ( ) 5 6
crack surface solid shell 1 2
(phantom nodes) coupling (BSM)

Coarse scale model
A
(c) A
Ghost atom Fine scale

(molecular statics)
Fine scale atom
Atom on the
crack surface
Continuum node
Gauss point
Fine scale model
Figure 7: a Schematic of a three-dimensional coupled continuum-atomistic model. b Mechanics of

coarse-scale domain modeled with solid shell element. c Fine-scale region showing the arrangement of
atoms in the diamond cubic lattice structure of Silicon. Pictures reproduced with permission from70.
procedure to couple a mesoscale discrete parti- material and size variables on the impact of
cle model and a macroscale continuum model design changes at the macroscopic level to the
to study the incompressible fluid flow, is pro- responses at the atomistic level.
posed in222. Li et al.223 developed a similar mul- However, extending the BSM to study
tiscale model, by combining the discrete element dynamic crack growth by adaptively adjusting the
method (DEM) at micro scale and Cosserat fine-scale domain and simultaneously avoiding
continuum modeling using the finite element spurious wave reflections at the ‘artificial’ inter-
method at macro scale, to simulate dynami- face is problematic. Budarapu et al.12 have devel-
cal responses in geo-structures composed of oped an adaptive multiscale method (AMM)
granular materials. A bridging scale method is to concurrently couple the atomistic domain
reported in224, for the analysis of localization with continuum by enhancing the bridging scale
problems. The micropolar-continuum model is method, to study crack propagation. They mod-
used to describe the localized deformation in a eled the coarse region based on the VAC model
small number of localized regions. A mathemati- and employed PNM to model crack in the coarse
cal framework of the bridging scale method and region. Furthermore, a meshless adaptive multi-
the time history kernel technique to impose the scale method for fracture (MAMMF) has been
dynamic interfacial boundary conditions are dis- reported in42, by coupling the atomistic domain
cussed in225. Farrell et al.35 employed the bridg- with DRKPM-based meshless method in the
ing scale method (BSM) to study intersonic crack continuum. Due to the absence of mesh in the
propagation, including the formation of a daugh- continuum, complex fracture patterns can be
ter cracks and the sudden acceleration of the captured with ease in the MAMMF. Recently, a
crack to a velocity exceeding the material shear solid-shell based three-dimensional concurrently
wave speed. They also proposed the non-reflect- coupled adaptive multiscale method has been
ing boundary conditions which can adequately introduced in70, to investigate the crack growth in
dissipate the strongly localized wave formed by thin-walled structures.
the Mach cone after the crack accelerates beyond
the material shear wave speed. The implementa-
4.2 Bridging domain method
tion algorithms as well as the development of a
Belytschko and Xiao6,7 have introduced the bridg-
time history kernel (THK) for the non-reflective
ing domain method for coupling the molecular
interface are discussed. A BSM-based model cou-
mechanics (molecular models at zero tempera-
pling the space–time Finite Element Method with
ture) and continuum models based on a domain
MD is developed in124 to simulate dynamic crack
decomposition technique. Coupling in the BDM
growth. A continuum-based sensitivity analysis
happens over a definite region (see Fig. 1c),
of two-dimensional continuum-atomistic mod-
known as ‘handshake’ domain. Unlike the BSM,
els using the bridging scale method is performed
in the BDM, the continuum region does not exist
in226. The authors correlated the influence of
352
(a) (b)
Figure 8: Schematic showing a close up of the region along the coupling boundary. Pictures reproduced
with permission from70.
in the atomistic domain. Second, in the BDM, domain in the bridging region B. In other words,
coupling between the atomistic and continuum in the handshake region the fine-scale displace-
regions is based on a linear energy weighting in ments must conform coarse-scale displacements.
the bridging domain and is enforced by Lagrange Therefore, considering the Lagrange multipliers,
multipliers. Therefore, the total energy of the the total Hamiltonian can be written as
system is a weighted contributions of the fine
and coarse models in the bridging domain. This HL = H + T g, (44)
is achieved through a scalar weight function, w where is Lagrange multipliers vector and g is
shown in Fig. 9, which is defined as unity out- the gap vector between coarse- and fine-scale
side, zero inside, the fine-scale domain and varies displacement. Derivation of explicit equations
smoothly in the blending region6: of motion for specific problems are explained in
several articles, see6,7,10,11,36,39,43,227,228.
∀X ∈ C \A

1
Belytschko et al.229 extended the BDM to cou-
w = [0, 1] ∀X ∈ B , (41) ple continua with molecular dynamics, where
0 ∀X ∈ A \C

the authors demonstrated the robustness of the
where C , B, and A correspond to the con- methodology by avoiding spurious wave reflec-
tinuum, bridging, and fine-scale domains, respec- tions at the molecular–continuum interface. A
tively. At any point X in the bridging domain, concurrent multiscale approach based on mul-
w can be computed by a normalized distance tigrid principles intended to solve large molecu-
function: lar dynamics systems is attempted in230. The
l(X)
authors estimated the effective stiffness matrix
w= , (42) of the coarse model by variational restriction
l0
of the effective stiffness matrix of the atomistic
where l(X) is the orthogonal projection of X on model. The influence of the time step and the dis-
the interior boundary of the coarse domain C cretization of Lagrange multipliers on spurious
and l0 is the length of this orthogonal projec- wave reflection are investigated in231 by includ-
tion to the boundary of the fine scale A, refer to ing a damping term in the fine-scale equations
Fig. 9. of motion. Guidault et al.36,227 have enhanced the
The governing equations of the coupled BDM by also enforcing the strain compatibility
model are derived from the Hamiltonian of the between the “atomistic” and continuum domains
coupled system, H, which is the sum of the Ham- in the bridging domain, which can be useful for
iltonians of each sub-domain: the development of error estimators to drive the
adaptive refinement of the coarse scale. However,
H = (1 − w)H A + wH C , (43)
little gain in accuracy is reported compared to the
where HAand HC
are the Hamiltonians of the much simpler L2 coupling. Xu et al.232 studied
fine- and coarse sub-domains, estimated as the the influence of the enforced constraints through
total potential of coarse- and fine-scale domains. Lagrange multipliers by modeling (1) exact non-
In the BDM, Lagrange multiplier method is diagonal Lagrange multiplier equations and
employed to constrain the coarse- and fine-scale (2) a diagonalized constraint form. Note that
353
the consistent constraint form conserves linear

w=0 w=1
momentum, angular momentum, and energy,
whereas the diagonalized constraint form dis-
l X
sipates energy. Therefore, the diagonalized form
is reported to be effective in suppressing spuri- l0
ous reflections at the interface. A variant of the
BDM for composite lattices through a ‘relaxed’
bridging domain method is discussed in233, where A
the atom set is divided into primary and second- C
ary atoms, and only the primary atoms are con-
strained in the coupling region. This will allow
Figure 9: The weight function in the handshake
the internal modes of the composite lattice to be domain in two dimensions.
relaxed, otherwise suppressed by the homogene-
ous continuum displacement field in the coupling
region.
(4) The zero KII criterion (Vanishing in-plane
Gracie et al.10 have extended the BDM for the
SIF (KII ) in shear mode for infinitesimally small
modeling of dislocations and cracks based on a
crack extension)242.
multiscale atomistic/continuum models. Fur-
In the multiscale methods mechanics of
thermore, they extended the multiscale model
crack growth around the crack tip are captured
to develop an adaptive concurrent multiscale
in the fine-scale region. Considering the atomis-
method228 for the dynamic simulation of dislo-
tic based multiscale methods, the crack growth
cations. Anciaux et al.234 have developed a BDM-
is identified based on the distance between the
based multiscale method to study the contact area
atoms. The orientation and branching are arrived
evolution of rough surfaces under normal load-
by following the path of the atoms on the crack
ing which can lead to the emergence of a strong
surface. This information is passed back to the
temperature gradient in the bridging zone. A gen-
coarse scale. Therefore, estimating the “length”
eralized bridging domain method is introduced
and orientation of crack growth and branching/
in235, based on independent weight functions to
joining cracks is more accurate and relatively easy
weight the material properties in the coarse- and
in multiscale methods. A two-scale approach to
fine-scale regions, followed by the force equilib-
simulate degradation and failure in polycrystal-
rium through imposing compensation forces esti-
line materials is proposed in243, where the macro-
mated by force and displacement compatibility
continuum is modeled using a three-dimensional
requirements. They tested the methodology on a
boundary element formulation in which the pres-
coupled continuum and discrete elements model.
ence of damage is formulated through an initial
stress approach to account for the local soften-
5 Fracture criteria and coarse graining ing in the neighborhood of points experiencing
A fracture criterion should determine whether degradation at the micro-scale. The two scales
a crack propagates/nucleates. It should further- are coupled by transferring the macro-strains to
more provide the orientation and “length” of the the micro-scale as periodic boundary conditions,
crack advancement, apart from whether or not while overall macro-stresses are obtained as vol-
cracks branch or join. Considering the non-line- ume averages of the micro-stress field.
arities and non-homogeneities around the crack To improve the computational efficiency, the
tip, estimating all the above details based on con- fine-scale region is adaptively adjusted as the
tinuum techniques alone is difficult, particularly crack propagates. The adaptivity scheme con-
when the “length” and direction of crack growth sists of adaptive refinement and coarse-graining
are not controlled66. This is because the fracture operations. In order to activate the adaptivity
criterion is often satisfied at several quadrature algorithm, the position of the crack tip in the
points in front of the crack tip, and reliable crite- fine-scale region (A) needs to be estimated. In
ria to estimate branching cracks are still missing. multiscale methods, when atomistic models are
Considering the approaches based on configura- used in the fine-scale region, the atoms on the
tional forces236–239, the four major cracking cri- crack surface are commonly identified based on
teria in LEFM are66 (1) Maximum hoop stress the energy criteria10 or centro symmetry param-
criterion or maximum principal stress criterion. eter (CSP)42,244,245. The crack tip is the ‘meeting
(2) Minimum strain energy density criterion240. point’ of atoms on either side of a crack sur-
(3) Maximum energy release rate criterion241. face12. Subsequently, crack growth is monitored
354
by comparing the location of the crack tip, in the of α. In an fcc lattice structure every atom α is
previous load step to the current load step. surrounded by 6 nearest neighbors (nnb). There-
fore, the CSP of the atom α in the fcc lattice is
given by
5.1 Energy criteria
The total potential energy of an atomistic system 3

is estimated as the sum of all bond potentials φα . CSPα = |rαβ + rα(β+3) |2 (47)
According to Eq. (4), the bond potential of a par- β=1
ticular atom depends on the distance between the
From Eq. (47), the CSP of an atom α in the fcc
atom (α) and its neighbors (β)10. In the initial
lattice is the summation of square of the total dis-
configuration, all the atoms are assumed to pos-
tance between the opposing neighbors. In other
sess the same potential energy. The initial crack is
words, the CSP of an atom in a periodic perfect
created by deleting the bonds between the atoms
lattice structure with symmetric atomic arrange-
and updating the neighbor list accordingly. Con-
ment is zero and the CSP values of the atoms on
tinuous increase in the external load leads to
the defect surface/stacking fault is not equal to
the stretching of the bonds of the atoms around
zero. This criterion is used to separate the atoms
the crack tip. Increase in bond length/distance
on the crack surface. Normalized CSP values for
between the atoms leads to increase in the system
various defects are listed in Table 1. From Table 1,
potential energy. A bond breaks when the bond
atoms on the crack surface can be distinguished
length reaches a certain threshold, transferring
as the atoms possessing normalized CSP values
the load to the immediate neighbors, see97 for
greater than or equal to 1.6881.
further details of bond elongation and rotation in
an initially notched Graphene system.
Therefore, the atoms around the crack tip 5.3 Adaptivity
possess the highest energy in the entire lattice. To improve the computational efficiency, the fine-
This is in agreement with the continuum theory scale region is adaptively enlarged with the defect
as well, where the stress concentration is observed propagation and the region behind the core of the
around the crack tip. Hence, potential energy defect (e.g., crack tip) is coarse grained. An adap-
provides an indication of the location of the tive concurrent multiscale methodology has been
crack tip. The energy criterion has been success- introduced in40 to handle the situations in which
fully applied to detect the locations of the crack both macroscopic and microscopic deformation
tip12 and the core of the dislocation228. Let EnHE be fields strongly interact near the tip of a crack.
the set of elements containing at least one atom The method is based on the balance between
with high potential energy, i.e. numerical and homogenization error; while the
first type of error states that elements should be
EnHE = {e ∈ EnA | energy of an atom in e > tolE }, refined in regions of high deformation gradients,
(45) the second implies that element size may not be
A E
where En is the set of total atoms and tol is the smaller than a threshold determined by the size
specified energy tolerance. As a guideline, tolE can of the unit cell representing the material’s micro-
be specified in the range of 15 and 30% higher structure. Adaptive refinement algorithms for 2D
than the energy of an atom in equilibrium in a peridynamic grids enhancing the concurrent mul-
perfect lattice. tiscale methods are discussed in246. They applied
the adaptive grid refinement to study dynamic
crack propagation in two-dimensional brittle
5.2 Centro symmetry parameter (CSP)
materials. An adaptive multiscale method cou-
The centro symmetry parameter of an atom α is
pling the space–time finite element method with
defined as follows76:
molecular dynamics is developed for the simula-
nb /2
n tion of dynamic fracture problems in124. Coupling
CSPα = |rαβ + rα(β+nnb /2) |2 , (46) between the fine- and coarse-scale simulation is
β=1 achieved with the introduction of a projection
operator and bridging scale treatment. An adap-
where rαβ and rα(β+nnb /2) are the distance between
tive multiscale methodology based on the Hill-
the atoms α and β and α and (β + nnb /2), respec-
Mandellemmainan FE2 sense is proposed to deal
tively, and nnb are the total number of nearest
with localized deformations in247. The displace-
neighbors of atom α. Consider an atom α in the
ment field of the fine-scale model was decom-
fine-scale region containing face-centered cubic
posed into a homogeneous part, fluctuations,
(fcc) lattice structure. Let β denote the neighbors
355
and a cracking part based on additional degrees

of freedom the crack opening in normal and tan- words, the minimum initial fine-scale region sat-
gential directions. Adaptive mesh refinement and isfying the above conditions is embedded within
coarsening schemes are proposed in248 for efficient a 3×3 discretization. Further details and applica-
computational simulation of dynamic cohesive tions to two- and three-dimensional crack growth
fracture. The adaptive mesh refinement consists problems are explained in12,42,70,245. The adaptiv-
of a sequence of edge-split operators, whereas the ity scheme consists of an adaptive refinement and
adaptive mesh coarsening is based on a sequence coarse-graining operations, as mentioned below:
of vertex-removal (or edge-collapse) operators.
The initial size of the fine-scale domain is cho- 1. Estimate the region in the coarse-scale
sen such that all the mechanics of crack growth domain C to be refined. A refinement
particularly around the crack tip are captured. operation involves the expansion of the fine-
Therefore, the initial domain size should be suf-
scale region by converting the estimated
ficient to surround the region ahead and behind
the crack tip. Furthermore, a large initial domain coarse region into a fine region, Fig. 10.
can lead to higher computational costs and the 2. Estimate the region in the fine-scale domain
crack tip may jump out of very small fine-scale A to be coarsened. In a coarse-graining
domains. Some of parameters that can influ- operation the coarse region is expanded by
ence size of the fine-scale region are as follows: converting the estimated fine region into a
(1) type of problem (static/dynamic), geometry coarse region, see Fig. 11.
and boundary conditions and (2) type and rate
of loading and the type of fracture (brittle/duc- In the above steps, when the sizes of the regions
tile). As a rule of thumb, a square domain is rec- refined and coarse grained are similar, the net
ommended: with areas behind and ahead of the change in the size of the fine-scale domain is
crack tip in the range of 20–25% and 80–75% of almost zero. As a result, the fine-scale region is
the total domain size, respectively, which leads to adaptively moved with the propagation of the
≈25% of the area behind the crack tip. defect.
Consider a fine-scale domain embedded
within the ’boundaries’ of the nodes/particles 5.3.1 Adaptive refinement
around the crack tip. The refinement algorithm The major steps of refinement (Fig. 10) proce-
should be activated sufficiently often such that a dure are listed for a multiscale method based on
buffer layer of elements/’regions’ is always main- an atomistic fine-scale model:
tained between the crack tip and the coupling
boundary. The ‘regions’ refer to the area/volume
1. Identify the region to be refined (ref ).
generated by connecting the immediate neigh-
2. Create and initialize the atoms in ref .
boring particles in meshless methods, such that
3. Identify and update the newly cracked
they resemble the elements in the mesh-based
atoms.
techniques. Second, to ensure that the refine-
4. Update the fine and coarse-scale regions.
ment operation is not activated in the first load
step itself, at least one layer of elements/regions is
Figure 10a shows the region identified for a
considered between the crack tip and the buffer
refinement operation. The fine-scale region after
element layer. Finally, the crack tip element layer
the refinement is depicted in Fig. 10b. Let the
is sandwiched by at least one layer of elements/
nodes/particles (before a refinement operation)
regions in the transverse direction. In other
in the fine, coarse, and completely cracked regions
split
be indicated by PnA, PnC, and Pn , respectively.
The region containing split elements indicates the
Table 1: Range of centro symmetry parameter completely cracked region. The steps of a refine-
for various defects, normalized by square of the ment operation are summarized as follows:
lattice parameter a02.
• Calculate the atoms on the crack surface based
Defect cspα /a02 Range �cspα /a02
on the CSP and store the regions containing
Perfect lattice 0.0000 cspα < 0.1 the atoms on the crack surface into the set
csp
Partial dislocation 0.1423 0.01 ≤ cspα < 2 Pn .
• Estimate the neighbours of the regions con-
Stacking fault 0.4966 0.2 ≤ cspα < 1 csp
taining the atoms on the crack surface in Pn
Surface atom 1.6881 cspα > 1 minA
and store them in Pn+1 .
356
• Calculate the regions to be refined, Pn+1

refine by • Identify the newly cracked particles in the
A
removing the fine-scale region Pn from the fine-scale region and update the split and tip
minA.
set Pn+1 nodes and the nodal connectivity table.
• Flag the regions to be refined and increase the
atomistic domain by creating the atoms in the A detailed algorithm of selecting the particles to
flagged elements. be refined, initializing the newly created atoms in
• Initialize the positions of the newly created the region identified for refinement and propa-
atoms through interpolation based on the gating the crack in the coarse-scale region in a
coarse-scale solution. multiscale framework, is explained in12,42.
• Update the fine and coarse regions after a
refinement operation. Update the neighbor
5.3.2 Adaptive coarse graining
list (nlistn+1 ) of the fine-scale atoms in the
The major steps for the coarse-graining operation
current load step (n + 1).
(Fig. 11) are as follows:
Area to be refined Crack
(a) (b)
Figure 10: Sketch of the adaptive refinement operation. a Flagged particles to be refined are hashed. b
Increased atomistic region after the refinement operation. Picture reproduced with permission from42.
Area to be coarsened Crack
(a) (b)
Figure 11: Schematic of the adaptive coarsening operation. a Flagged particles to be coarsened are
hashed. b Reduced atomistic region after the coarsening operation. Picture reproduced with permission
from42.
357
1. Identify the fine-scale region to be coarse 2. Identify the regions containing atoms on the
grained (coa). crack surface, based on the positions of the
2. Estimate the equivalent coarse-scale region atoms on the crack surface and the positions
of coa, refer Sect. 5.4. of the particles/nodes of the background
3. Delete the atoms in the region to be coars- discretization, see Fig. 12b.
ened. 3. Estimate the normal and center of gravity
4. Update the fine- and coarse-scale particles/ (CoG) of the atoms on the crack surface.
nodes. Calculate the effective CoG of a crack region
by averaging the CoGs of the atoms on the
crack surface in the considered crack region.
The process of an adaptive coarse-graining 4. Approximate the crack path in each crack
operation is explained in Fig. 11. Let PnCS be the region by joining the effective normal and
regions containing atoms on the crack surface at CoG of the atoms on the crack surface, refer
load step n. Let PnBA be the regions lying in the Fig. 12d and Sect. 5.4.1.
fine-scale domain and attached to the coupling 5. Estimate the nodes or particles on either side
‘boundary’. The particles/nodes to be coars- of the crack surface or around the tip, see
ened are the particles/nodes which are in both Fig. 12c.
set PnCS and the set PnBA in front of the crack tip,
Pncoarsen = PnCS ∩ PnBA. The steps of an adaptive
coarsening operation are as follows: 5.4.1 Crack surface orientation
Consider a deformed configuration of the fine-
• Estimate and store the regions containing the scale model, superimposed with a discretized
elements on the crack surface (far away from coarse-scale model as shown in Fig. 12a. The
the crack tip) into PnLE. atoms in the fine region can be separated into
• Find the fine-scale regions attached to the small rectangular cells surrounded by four nodes/
coupling boundary, PnBA. particles in the coarse region. The center of grav-
• The regions to be coarse grained (Pn+1
coarsen) are
ity of a cell containing the atoms on the crack
LE
given by Pn ∩ Pn . BA
surface can be calculated by averaging the posi-
• Flag the regions to be coarse grained and tions of center of gravities of the atoms on the
delete the atoms in the flagged regions. cog
crack surface (rα ) in that cell245:
• Update the particle/nodal set in the fine- and ncacr cog
coarse-scale regions and the neighbor list of cog rα
rcell = α=1cs , (48)
the fine-scale atoms, after a coarsening opera- n
tion. cog
where rcell is the approximated position of the
center of gravity of the atoms on the crack sur-
More details can be found in references12,42. face and ncs are the total number of atoms on
the crack surface, in a crack region. The normal
5.4 Efficient coarse‑graining techniques of the approximated crack surface in the crack
Upscaling the fracture-related material informa- region is computed as the average of the normals
tion from the fine scale to the coarse scale is a of the atoms on the crack surface:
major difficulty in multiscale methods for frac- ncs cog
cog nα
ture, particularly for complex crack patterns. ncell = α=1cs , (49)
Belytschko et al.249 developed a coarse-graining n
cog
approach named multiscale aggregating dis- where ncell is the normal vector of the approxi-
continuity method. A robust and simple coarse- mated crack surface in a crack region. There-
graining technique in the context of multiscale fore, the crack surfaces in the crack regions are
modeling for fracture is developed by Budarapu obtained based on the planes passing through
cog
et al.245 The major steps in245 to develop an equiv- rcell , whose normals are estimated from Eq. (49).
alent model of the A Finally, the approximated crack surface in the CG
def , the coarse-graining (CG)
method (Fig. 12), are as follows: model are obtained by joining the crack surfaces
in each crack region.
In order to generate a smooth and continu-
1. Determine the atoms on the crack surface,
ous crack surface in the CG domain, the start/
e.g., using the CSP. end positions of the crack surfaces on the verti-
cal edges of the crack regions are averaged, as
358
illustrated in the schematic Fig. 13. As a rule of 6 Computer implementation

thumb, a cell containing at least 12 atoms on In this section, we discuss the numerical imple-
the crack surface is observed to be considered mentation details of a three-dimensional multi-
as crack region245. Therefore, the minimum size scale method for fracture. The atomistic model
of the cell can be adopted as 13 times the lattice in the fine region is assumed to be modeled using
parameter. The cell size or the size of fine-scale LAMMPS. The whole computational framework
domain in general could be determined by a-pos- is developed in MATLAB, where the LAMMPS is
teriori error estimators. An example of generation triggered through system command.
of a continuous crack surface in the coarse region
is demonstrated in Fig. 13. Consider the vertical
6.1 Codes and algorithms
edge containing points C, D, E, and F. The points
Consider a coarse-scale model implemented in
D and E correspond to end points of two crack
MATLAB coupled to a fine-scale model in the
surfaces and the points C and F are the starting
LAMMPS software. Since the LAMMPS software
points of new crack surfaces. The largest distance
can be triggered from MATLAB, a versatile and
between these points is the distance between the
robust multiscale strategy can be developed in the
points C and F which is larger than the domain of
MATLAB frame work. To develop such numerical
influence. Thus there exists more than one point
methodology, the system command in MAT-
on the equivalent crack surface on this particular
LAB, which triggers an executing system opera-
edge. The total number of points on the equiva-
tion, is used as described below70:
lent crack surface on this vertical edge can be
estimated by recursively checking if the distance system(‘lmp_mpi − in . . . / . . . /input_file_name − loglog.ini′ );
between the neighbors of points C, D, E, and F (50)
falls within the domain of influence. Figure 14
where ‘lmp_mpi’ indicates the LAMMPS execut-
shows the equivalent coarse-grained model of
able file generated by compiling the parallel ver-
an atomistic model245. The deformed configura-
sion of the LAMMPS code. The command ‘.../.../
tion of the atomistic model for a dynamic double
input_file_name’ is used to identify the exact
edge crack propagation after 108 pico-seconds is
location of the input file.
shown in Fig. 14a. The corresponding equivalent
coarse-grained model is shown in Fig. 14b.
(c)
(a) (b) (d)
Figure 12: Schematic of meshless equivalent coarse-scale model. a Meshfree particles superimposed

on the atomistic model, b regions containing the atoms on the crack surface are highlighted along with
the normals of the crack surface in each region, c calculation of level sets and d approximation of crack
surface by joining the crack path in the regions containing the crack. Picture reproduced from245 with per-
mission.
359
6.1.1 Implementation in LAMMPS the current step is represented by ‘ubary’(=0.05)

Algorithm 1 describes the key steps in a variable. The atoms in the groups ‘top’ and ‘bot’
LAMMPS input file to estimate the atom posi- are uniformly displaced by an amount of ‘ubary’
tions through energy minimization. LAMMPS through the ‘displace_atoms’ command. Com-
commands are highlighted in blue. The load is mand ‘fix’ helps in maintaining the bound-
assumed to be prescribed in several steps through ary conditions of the specified group of atoms
‘nsteps’(=100) variable, and the atom positions through ‘setforce’ option, where a ‘NULL’ value
at each load step are estimated based on the mini- indicates no constraint and ‘0.0’ denote a con-
mum energy. ‘step_c’ indicates the current load straint in that direction. For example, the line “fix
step and the amplitude of the displacement in 2 top setforce NULL 0.0 NULL”, reads
Figure 13: Schematic of averaging the approximated individual crack surface orientation in each crack
region, to generate a smooth continuous equivalent crack surface. Picture reproduced from245 with per-
mission.
800
700
Y − coordinate (angstroms)
600
500
400
300
200
100
−100
0 200 400 600 800 1000 1200
X − coordinate (angstroms)
(a) (b)
Figure 14: Development of an equivalent coarse-scale model of a given fine-scale model, for a dynamic
crack propagation of double edge crack model. a Deformed configuration after 108 pico-seconds along
with the highlighted atoms on the crack surface, crack regions, and their normals. Approximated crack
surface showing the corresponding approximated equivalent crack surfaces. Pictures reproduced from245
with permission.
360
as follows: fix number 2, where the ‘top’ group of Fig. 16b, indicates a corresponding drop in the
atoms are constrained along the ‘y’ direction and force. After breaking the initial bonds, the mate-
free to move in the ‘x’ and ‘z’ directions. Finally, rial resisting the external load will be continuously
the energy minimization is carried out through diminishing as the material separation continues.
the ‘minimize’ command until the convergence
of atom positions within the prescribed limits.
6.1.2 Multiscale model
The process is repeated in the next load step. The
The major steps of multiscale model for fracture
results are stored in the specified file through the
are summarized in Algorithm 2. The split and tip
‘dump’ command before proceeding to the next
elements are identified in the discretized domain
load step.
based on the geometry of the crack(s). This is
Consider the simulation of punching a hole in followed by generating the initial configuration of
a rectangular panel through molecular dynamics. the atomistic model including initial notches (if
This is achieved by specifying a uniform out-of- applicable). In the third step, a for loop applies
plane displacement along the z-direction to the the boundary conditions on the coarse-scale in
upper side of the plate in a specific area. Due to several load steps providing the coarse-scale solu-
symmetry, a quarter of the plate is considered. tion uIC, in each step. The ghost atom positions
The hole is punched using a quarter circle with are interpolated from the coarse-scale solution
a prescribed radius. The quarter circle is further using Eq. (40). They are the boundary conditions
extruded to a quarter cylinder along the thick- for the fine scale. The LAMMPS executable can
ness direction. Atoms on the upper side of the now be triggered again to minimize the poten-
quarter cylinder portion are subjected to uniform tial energy of atomistic domain by fixing the
displacements along the z-direction. Figure 15a–c updated ghost atom positions. In each load step,
shows the initial and deformed configuration the latest atom positions at the end of the energy
after 100 and 1530 load steps, respectively. The minimization along with their energy and centro
entire material is separated from the plate mate- symmetry parameter (CSP) are stored an output
rial after 1530 steps. The evolution of the potential file. The crack tip is identified by using either the
energy versus the strain is plotted in Fig. 16a. The energy or the CSP criteria. The adaptivity scheme
potential energy fluctuates after exceeding a strain is activated if the crack tip location is close to the
of 0.48. The load–displacement curve plotted in boundary of the atomistic domain.
361
Figure 15: Deformed configurations based on the atomistic model, during punching a hole in Silicon.
Distribution of the potential energy at various instances of the punching process. a Initial, b after 100 load
steps, and c after 1530 steps. Pictures reproduced with permission from70.
× 10 5
-3.36
-3.365
Potential energy (eV)
-3.37
-3.375
-3.38
-3.385
-3.39
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
Strain along the z direction ( zz
)
(a) (b)
Figure 16: a Distribution of the potential energy with strain and b load-displacement diagram generated
based on the MD simulations of punching a hole in a rectangular panel. Pictures reproduced with permis-
sion from70.
Consider the simulation of Mode I crack simulation is shown in Fig. 17d, where almost a
growth of two through-the-thickness edge cracks, complete merging of the two cracks and hence
located in the middle of the plate. The cracks in the separation of the fine-scale region into two
the coarse-scale domain are modeled by employ- parts can be noticed. Some key contributions
ing the phantom node method. Displacement towards the multiscale methods for fracture in
loads are prescribed on the top and bottom row the past two decades are summarized in Table 2.
of atoms. Deformed configurations after 28 and
38 load steps are plotted in Fig. 17a, b, respec-
7 Future prospects and conclusions
tively. Cracks propagate in the opposite directions
Most of the problems in real-time involve mul-
after 28 load steps. Figure 17c shows the cou-
tiple field and disparate time and length scales.
pled model after an adaptive refinement after 39
Based on the significant advancement of the
load steps. Since the available space between the
multiscale methods in the past two decades, mul-
initial fine-scale regions in Fig. 17a is small, the
tiphysics multiscale methods are rapidly grow-
two fine-scale regions are merged in the adaptive
ing to simulate fracture in various applications,
refinement. Simultaneously, the adaptive coars-
such as heterogeneous porous medis and/or
ening scheme coarse grains the fine-scale regions
hydralic fracture253, polycrystalline254 and com-
behind the crack tips. The deformed configu-
posite materials255,256 design, batteries257,258, nano
ration of the multiscale model at the end of the
362
Figure 17: Crack propagation of a double-edge notched specimen. a Deformed configuration after 28

load steps. Activation of adaptive refinement and coarse graining of the fine-scale region as the crack
grows. b Deformed configuration after 38 load steps, before refinement and c after an adaptive refine-
ment operation after 39 load steps. Adaptive refinement and coarse-graining algorithms (see12) are acti-
vated after 39 load steps as the cracks grow. As a result, the two fine-scale regions are merged after
39 load steps and the combined fine-scale region is adaptively coarse grained, as shown in d. Pictures
reproduced with permission from70.
materials259, bio applications260, to name a few. 7.1 Quantum mechanics and molecular

In this section, we outlined the current status of mechanics
latest multiscale techniques based on quantum Estimation of mechanical and fracture prop-
mechanics, peridynamics, and techniques for bio- erties of the nano scale structures such as car-
logical applications, apart from some comments bon nanotubes, two-dimensional materials like
on future prospects. Graphene and MoS2 through experiments is
363
extremely challenging261. Hence, as an alternative approach by concurrently combining the QM,

to experiments, quantum mechanics (QM)based MM, and CG techniques. They identified two
techniques can be primarily employed to predict distinct interfaces, MM/CG and QM/MM, where
the fracture, bond breaking, and bond formation the QM domain does not interact with CG
in particular262, and when using the simulation domain and MM and CG regions are coupled by
techniques that consider electrons, i.e., quantum quasi-continuum.
modeling is indispensable. In spite of significant
progress in the computational resources, the
7.2 Peridynamics based multiscale
computational expenses of simulations based on
methods
QM still remains a challenging task. Therefore,
Peridynamics266–269 is a nonlocal computational
MD simulations are computationally more eco-
formulation of continuum mechanics, equivalent
nomical compared to that of QM. However, due
to a coarse-grain model in multiscale perspec-
to small dimensions of atoms, full-scale atomis-
tive. The main difference between the continuum
tic models for day-to-day engineering calcula-
mechanics and the peridynamics is the nonlocal
tions are prohibitively expensive. In this context,
interaction between material points. Consider
multiscale methods coupling quantum mechan-
a fixed material point rα in the current con-
ics, molecular mechanics (MM), coarse graining
figuration, which can interact with neighbor-
(CG), and continuum mechanics (CM) tech-
ing particles rβ within a compact support called
niques seem to be efficient, coupling the advan-
as horizon, which is similar to the concept of
tages and surpassing the disadvantages of each
‘domain of influence’ in the meshless methods
technique individually.
or ‘cutoff range’ in molecular dynamics. The
ONIOM263 is a technique coupling QM and
nonlocal equilibrium equations of motion can
MM based on overlapping domain scheme.
be derived by considering the balance of linear
The basic idea of ONIOM scheme is to apply
momentum at the material point rα, as men-
the coarse-scale model like MM to the entire
tioned below180:
domain and the more accurate QM-based fine-
scale model in the critical regions where bond ρα r̈α = L(Rα , t) + ρα b(Rα ), (52)
breaking is expected. A higher order correction,
where ρα is the average density of the material at
estimated as the difference between the QM and
point α
MM energies of the fragment domain, can be
applied to cancel the effect of the complete MM

L(Rα , t) = (Tα �Rβ − Rα � − Tβ �Rα − Rβ �)dVβ
model assumption in the coupling domain that Hα
is overlapped by the QM model. Therefore, the
(53)
total energy of the system can be written263 as
is the nonlocal stress divergence vector acting on
follows:
the α th material point by neighboring macroscale
QM points β, equivalent to a local divergence term in
E(xα , cα ) = E MM (xα ) + EF (xα , cα ) − EFMM (xα ), continuum mechanics180.
(51) A peridynamics based multiscale method has
been developed in180, using the multiscale micro-
where E MM is the energy of the MM domain,
QM morphic molecular dynamics (MMMD) theory
EF is the energy of the fine-scale sub-domain
to couple the molecular dynamics in the fine-
calculated based on QM, EFMM is the energy
scale region. To address the issue of wave reflec-
estimated using MM in the fine-scale region, xα
tion on the interface, the authors proposed a
is the vector of current atom position, and cα is
filter by turning on and off the MMMD dynamic
the set of basis function coefficients used in the
equations at different scales. A coupled model
representation of the electronic wave function.
embedding peridynamics within a molecular
On the other hand, the quantum to molecular
dynamics code is available at270. A peridynamics
mechanical overlapping domain (QtMMOD)
based hierarchical multiscale modeling scheme
method264 requires a partial overlap between the
coupling peridynamics with the atomistic model
MM and the QM sub-domains. Therefore, the
has been employed in271, to model a complex het-
QM model is used in the interesting area, whereas
erogeneous polymer, ultra high molecular weight
the underlying MM model is employed on the
polyethylene (UHMWPE). Refer to179,252 for a
sub-domain excluding the interesting region.
multiscale strategy coupling peridynamics with
In other words, the MM domain does not exist
continuum based finite element method. Differ-
in the entire domain as in the case of ONIOM
ent methods are compared in272, to estimate the
method. Park et al.265 propose another multiscale
364
Table 2: Key contributions on multiscale methods for fracture based on BSM and BDM techniques.
No. Address Title Comments
1 Wagner et al., J Comp Phy, Coupling of atomistic and contin- Introduction to BSM
190:249–279, 2003 uum simulations using a bridging
scale decomposition32
2 Park et al., Phil Mag, 85:79–113, The Bridging Scale for Two- Wave relections and time history
2005 Dimensional Atomistic/Continuum kernel in BSM
Coupling33
3 Tang et al., IJNME, 65:1688–1713, A mathematical framework of the Modified interfacial conditions
2006 bridging scale method225 based on THK
4 Farrell et al., IJNME, 71:583–605, Implementation Aspects of the Computer implementation
2007 Bridging Scale Method and algorithms and shear dominant
Application to Intersonic Crack failure using BSM.
Propagation35
5 Budarapu et al., CMAME, 319:338– Concurrently coupled solid shell- Solid shell-based multiscale
365, 2017 based adaptive multiscale method method for adaptive crack
for fracture70 growth, using BSM.
6 Belytschko et al., CMAME, A bridging domain method for BDM to simulate dynamic frac-
193:1645–1669, 2004 coupling continua with molecular ture.
dynamics229
7 Gracie et al., IJNME, 86:575–597, Adaptive Continuum-Atomistic Adaptive XBDM for dislocation
2011 Simulations of Dislocation Dynam- dynamics.
ics228
8 Talebi et al., Comp Mech, 53:1047– A Computational Library for Multiscale framework for 3D
1071, 2014 Multiscale Modelling of Material dynamic fracture using XBDM.
Failure11
9 Miller et al., MSMSE, 17:053001, A unified framework and perfor- Comparison of accuracy and effi-
2009 mance benchmark of fourteen ciency of 14 multiscale methods
multiscale atomistic/continuum on a test problem.
coupling methods250
10 Nair et al., JMPS, 59:2476–2487, ACoupled quantum-continuum Coupled quantum-continuum
2011 analysis of crack tip processes in analysis of crack.
aluminum251
11 Liu et al., CMAME, 245:163–175, A coupling approach of discretized Coupled FEM and Peridynamics
2012 peridynamics with finite element model.
method252
12 Giovanardi et al., CMAME, 2017 A hybrid XFEM-Phase field (Xfield) Coupled XFEM-Phasefield method.
method for crack propagation in
brittle elastic materials218
tangent-stiffness matrices in a massively parallel prototypes they are mimicking. The main rea-
computational peridynamics code. son being the weak junctions between stiff and
compliant phases, which need to be optimized
according to the intended functions of the com-
7.3 Multiscale methods for biological
posite material273.
applications
Investigating biological system mechanics at
Many biological materials, such as nacre, tooth,
the smallest scale does not always provide a com-
and bone are composite materials made up of
plete picture274. Therefore, understanding the
stiff brittle ceramics and compliant organic
influence of multiphasic interfaces and hierarchi-
materials like polymer. Compared to their con-
cal organizations across length scales on macro-
stituents, natural organic/inorganic composites
scale properties in natural systems will help in
exhibit much enhanced strength and toughness
developing a materials-by-design approach for
properties. Based on this inspiration, several bio-
novel engineering materials, such as nanocom-
mimetic composites are proposed in an attempt
posites with tailored interfaces and programmed
to synthesize materials with superior mechanical
microstructures, e.g., the brick-and-mortar
properties. However, most current synthetic com-
arrangement of stiff filler and soft matrix phases
posites have not exhibited their full potential of
observed in nacre275. Niebel et al.276 conducted
property enhancement compared to the natural
365
experimental and numerical studies to under- interdisciplinary collaboration among differ-

stand the influence of the polymer properties on ent groups is required to accurately predict and
the mechanics of nacre-like composites contain- understand the physics/mechanics across the
ing an intermediate fraction of mineral phase and scales. This is a good sign, which helps not only
reported that the stiffer polymers can increase to quickly understand the fundamental mechan-
the strength of the composite by reducing stress ics of failure, but also to bring in revolutionary
concentrations in the inorganic scaffold. A finite changes in the material design and analysis.
element based analysis is carried out in273, to
estimate the improvement in the mechanical
Acknowledgements
properties of nacre like biomimetic composites.
PRB acknowledge the funding from the Euro-
Awaja et al.277 summarized the recent develop-
pean Research Council (ERC), Grant No.
ments on the topics of cracks and microcracks
306622 through the ERC Starting Grant “Multi-
initiation and propagation in polymer struc-
field and multi-scale Computational Approach
tures along with the techniques for detection
to Design and Durability of PhotoVoltaic
and observation. Moreover, repair of cracks and
Modules”—CA2PVM.
microcracks through bio-mimetic self-healing
techniques is also discussed along with surface Received: 20 January 2017 Accepted: 19 July 2017
active protection. Published online: 11 October 2017
Mechanisms of energy dissipation in struc-

tural molecules at nanoscales can be estimated
based on the sacrificial bonds and hidden length References
(SBHL). The presence of SBHL leads to greater 1. Tadmor EB, Ortiz M, Phillips R (1996) Quasicon-
fracture toughness as compared to the materials tinuum analysis of defects in solids. Philos Mag A
without such features. The increase in interface 73(6):1529–1563
toughness as a function of polymer density and 2. Abraham FF, Walkup R, Gao H, Duchaineau M,
number of sacrificial bonds has been investigated DeLaRubia TD, Seager M (2002) Simulating materials
in278, based on the mechanical properties of the failure by using up to one billion atoms and the world’s
polymeric system. Tessellation is a structural fastest computer: work-hardening. Proc Nat Acad Sci
motif involving periodic soft and hard elements 99(9):5777–5782
arranged in series which appears in a vast array of 3. Buehler MJ, Hartmaier A, Gao H, Duchaineau M,
invertebrate and vertebrate animal biomaterials. Abraham FF (2004) Atomic plasticity: description and
Tessellation of a hard, continuous surface, con- analysis of a one-billion atom simulation of ductile
nected by a softer phase results in maximization materials failure. Comput Methods Appl Mech Eng
of material toughness, with little expense to stiff- 193(48–51):5257–5282
ness or strength279. 4. Liu WK, Su H, Belytschko T, Li S, Chang CT (2000)
Dental enamel is a hybrid material consisting Multi-scale methods. Int J Numer Methods Eng
of brittle fibers and compliant organic materials 47:1343–1361
like protein matrix. Enamel exhibits high frac- 5. Miller RE, Tadmor EB (2002) The quasicontinuum
ture toughness and stiffness due to a complex method: overview, applications and current directions.
hierarchical and graded microstructure, opti- J Comput Aided Mater Des 9(3):203–239
mally organized from nano to macroscale. The 6. Belytschko T, Xiao SP (2003) Coupling methods for
deformation and damage behavior of the fibrous continuum model with molecular model. Int J Multi-
microstructure is studied in280 using a 3D RVE scale Comput Eng 1(1):115–126
and continuum damage mechanics model cou- 7. Xiao SP, Belytschko T (2004) A bridging domain
pled to hyperelasticity for modeling the initiation method for coupling continua with molecu-
and evolution of damage in the mineral fibers as lar dynamics. Comput Methods Appl Mech Eng
well as protein matrix. Ural et al.281 simulated the 193(17–20):1645–1669
bone fracture based a multiscale method using 8. Liu WK, Park HS, Qian D, Karpov EG, Kadowaki H,
cohesive finite elements. A failure mode transi- Wagner GJ (2006) Bridging scale methods for nano-
tion in nacre and bone-like materials has been mechanics and materials. Comput Methods Appl Mech
demonstrated in282. Eng 195(13–16):1407–1421
To summarize, based on the rapid progress 9. Belytschko T, Loehnert S, Song JH (2008) Multiscale
in multiscale methods, the expensive and time- aggregating discontinuities: a method for circumvent-
consuming experiments can be avoided in future ing loss of material stability. Int J Numer Methods Eng
either completely or partially. Moreover, since 73(6):869–894
the real-time problems involve multiphysics, an
366
10. Gracie R, Belytschko T (2008) Concurrently cou- fiber sic/ti composite materials. Comput Methods Appl
pled atomistic and XFEM models for dislocations and Mech Eng 183:309–330
cracks. Int J Numer Meth Eng 78(3):354–378 26. Feyel F, Chaboche JL (2001) Multi-scale non linear FE2
11. Talebi H, Silani M, Bordas SPA, Kerfriden P, Rabczuk T analysis of composite structures: damage and fiber size
(2014) A computational library for multiscale model- effects. In: Saanouni K (ed.) Numerical Modelling in
ling of material failure. Comput Mech 53(5):1047–1071 Damage Mechanics - NUMEDAM00. Revue Eur Elem
12. Budarapu PR, Gracie R, Bordas SPA, Rabczuk T (2014) Finis 10:449–472
An adaptive multiscale method for quasi-static crack 27. Feyel F, Chaboche JL (2003) A multilevel finite element
growth. Comput Mech 53(6):1129–1148 method (FE2) to describe the response of highly non-
13. Shenoy VB, Miller RE, Tadmor EB, Rodney D, Phillips linear structures using generalized continua. Comput
R, Ortiz M (1999) An adaptive finite element approach Methods Appl Mech Eng 192:3233–3244
to atomic-scale mechanics-the quasicontinuum 28. Talebi H, Zi G, Silani M, Samaniego E, Rabczuk T
method. J Mech Phys Solids 47(3):611–642 (2012) A simple circular cell method for multi-level
14. Beex LAA, Peerlings RHJ, Geers MGD (2011) A qua- finite element analysis. J Appl Math. ID: 526846
sicontinuum methodology for multiscale analysis of 29. Silani M, Ziaei-Rad S, Talebi H, Rabczuk T (2014) A semi-
discrete microstructural models. Int J Numer Methods concurrent multiscale approach for modeling damage in
Eng 87:701–718 nanocomposites. Theoret Appl Fract Mech 74:30–38
15. Sun Y, Peng Q, Lu G (2013) Quantum mechanical mod- 30. Zhu H, Wang Q, Zhuang X (2016) A nonlinear semi-
eling of hydrogen assisted cracking in aluminum. Phys concurrent multiscale method for fractures. Int J
Rev B 88:104109 Impact Eng 87:65–82
16. Beex LAA, Kerfriden P, Rabczuk T, Bordas SPA (2014) 31. Zhuang X, Wang Q, Zhu H (2017) Multiscale modelling
Quasicontinuum-based multiscale approaches for of hydro-mechanical couplings in quasi-brittle materi-
plate-like beam lattices experiencing in-plane and out- als. Int J Fract. doi:10.1007/s10704-016-0139-1
of-plane deformation. Comput Methods Appl Mech 32. Wagner GJ, Liu WK (2003) Coupling of atomistic and
Eng 279:348–378 continuum simulations using a bridging scale decom-
17. Beex LAA, Peerlings RHJ, Geers MGD (2014) A multi- position. J Comput Phys 190(1):249–279
scale quasicontinuum method for lattice models with 33. Park HS, Karpov EG, Liu WK, Klein PA (2005) The
bond failure and fiber sliding. Comput Methods Appl bridging scale for two-dimensional atomistic/contin-
Mech Eng 269:108–122 uum coupling. Philos Mag 85(1):79–113
18. Hajibeygi H, Karvounis D, Jenny P (2011) A hierarchi- 34. Dhia HB (2006) The Arlequin method: a partition of
cal fracture model for the iterative multiscale finite vol- models for concurrent multiscale analyses. In: Chal-
ume method. J Comput Phys 230:8729–8743 lenges in computational mechanics workshop, Cachan,
19. Paggi M, Wriggers P (2012) Stiffness and strength of France, 10–12 May, 2006
hierarchical polycrystalline materials with imperfect 35. Farrell DE, Park HS, Liu WK (2007) Implementation
interfaces. J Mech Phys Solids 60:557–572 aspects of the bridging scale method and application to
20. Paggi M, Corrado M, Rodriguez MA (2013) A multi- intersonic crack propagation. Int J Numer Methods Eng
physics and multi-scale numerical approach to microc- 71:583–605
racking and power-loss in photovoltaic modules. Com- 36. Guidault PA, Belytschko T (2009) Bridging domain
pos Struct 95:630–638 methods for coupled atomistic continuum mod-
21. Wudtke I, Talebi H, Silani M, Werner F (2015) A hierar- els with l 2 or h1 couplings. Int J Numer Methods Eng
chical multi-scale approach to mechanical characteriza- 77(4–5):1566–1592
tion of heat affected zone in welded connections. Com- 37. Broughton J, Abraham FF, Bernstein N, Kaxiras E
put Mater Sci 96:396–402 (1999) Concurrent coupling of length scales: methodol-
22. Lawrimore WB, Paliwal B, Chandler MQ, Johnson KL, ogy and application. Phys Rev B 60(4):2391–2403
Horstemeyer MF (2016) Hierarchical multiscale mod- 38. Tang S, Kopacz AM, O’Keefe SC, Olson GB, Liu WK
eling of polyvinyl alcohol/montmorillonite nanocom- (2013) Concurrent multiresolution finite element:
posites. Polymer 99:386–398 formulation and algorithmic aspects. Comput Mech
23. Lyu D, Fan H, Li S (2016) A hierarchical multiscale 52(6):1265–1279
cohesive zone model and simulation of dynamic frac- 39. Talebi H, Silani M, Bordas SPA, Kerfriden P, Rabczuk T
ture in metals. Eng Fract Mech 163:327–347 (2013) Molecular dynamics/XFEM coupling by a three
24. Feyel F (1999) Multiscale FE2 elastoviscoplastic dimensional extended bridging domain with applica-
analysis of composite structures. Comput Mater Sci tions to dynamic brittle fracture. Int J Multiscale Com-
16(1–4):344–354 put Eng 11(6):527–541
25. Feyel F, Chaboche JL (2000) FE2 multiscale approach 40. Vernerey FJ, Kabiri M (2014) Adaptive concurrent
for modeling the elastoviscoplastic behavior of long multiscale model for fracture and crack propagation
367
in heterogeneous media. Comput Methods Appl Mech 55. Sheng Y, Sousani M, Ingham D, Pourkashanian M
Eng 276:566–588 (2015) Recent developments in multiscale and mul-
41. Wu J, Zhang H, Zheng Y (2015) A concurrent multi- tiphase modelling of the hydraulic fracturing process.
scale method for simulation of crack propagation. Acta Math Probl Eng 2015:729672
Mech Solida Sin 28(3):235–251 56. Liu Y, Filonova V, Hu N, Yuan Z, Fish J, Yuan Z, Belytschko T
42. Yang S-W, Budarapu PR, Mahapatra DR, Bordas SPA, (2014) A regularized phenomenological multiscale dam-
Zi G, Rabczuk T (2015) A meshless adaptive multiscale age model. Int J Numer Methods Eng 99:867–887
method for fracture. Comput Mater Sci 96B:382–395 57. Paggi M, Wriggers P (2011) A nonlocal cohesive zone
43. Talebi H, Silani M, Rabczuk T (2015) Concurrent mul- model for finite thickness interfaces—Part I: math-
tiscale modeling of three dimensional crack and dislo- ematical formulation and validation with molecular
cation propagation. Adv Eng Softw 80:82–92 dynamics. Comput Mater Sci 50:1625–1633
44. Ghosh S, Lee K, Moorthy S (1994) Multiple scale analy- 58. Paggi M, Wriggers P (2011) A nonlocal cohesive zone
sis of heterogeneous elastic structures using homogeni- model for finite thickness interfaces—Part II: FE imple-
zation theory and voronoi cell finite element method. mentation and application to polycrystalline materials.
Int J Solids Struct 32:27–62 Comput Mater Sci 50(5):1634–1643
45. Kouznetsova V, Geers MGD, Brekelsmans WAM (2002) 59. Paggi M, Reinoso J (2015) An anisotropic large dis-
Multi-scale constitutive modeling of heterogeneous placement cohesive zone model for fibrillar and crazing
materials with a gradient-enhanced computational of interfaces. Int J Solids Struct 69–70:106–120
homogenization scheme. Int J Numer Methods Eng 60. Shojaei A, Li G, Fish J, Tan PJ (2014) Multi-scale con-
54:1235–1260 stitutive modeling of ceramic matrix composites by
46. Guidault PA, Allix O, Champaney L, Navarro JP (2007) continuum damage mechanics. Int J Solids Struct
A two-scale approach with homogenization for the 51:4068–4081
computation of cracked structures. Comput Struct 61. Greco F, Leonetti L, Luciano R, Blasi PN (2016) An
85:1360–1371 adaptive multiscale strategy for the damage analysis of
47. Özdemir I, Brekelmans WAM, Geers MGD (2008) masonry modeled as a composite material. Compos
FE2 computational homogenization for the thermo- Struct 153:972–988
mechanical analysis of heterogeneous solids. Comput 62. Nemat-Nasser S, Hori M (1993) Micromechanics:
Methods Appl Mech Eng 192:602–613 overall properties of heterogeneous materials. Elsevier,
48. Verhoosel CV, Remmers JJC, Gutiérrez MA, de Borst Amsterdam
R (2010) Computational homogenization for adhesive 63. Verhoosel CV, Remmers JJC, Gutiérrez MA (2010) A
and cohesive failure in quasi-brittle solids. Int J Numer partition of unity-based multiscale approach for mod-
Methods Eng 83(8–9):1155–1179 elling fracture in piezoelectric ceramics. Int J Numer
49. Nguyen VP, Lloberas-Valls O, Stroeven M, Sluys Methods Eng 82:966–994
LJ (2011) Homogenization-based multiscale 64. Oliver J, Caicedo M, Huespe AE, Hernández JA, Roubin
crack modelling: From micro-diffusive damage to E (2017) Reduced order modeling strategies for com-
macro-cracks. Comput Methods Appl Mech Eng putational multiscale fracture. Comput Methods Appl
200(9–12):1220–1236 Mech Eng 313:560–595
50. Nguyen VP, Lloberas-Valls O, Stroeven M, Sluys LJ 65. Xu M, Gracie R, Belytschko T (2009) Bridging the
(2012) Computational homogenization for multiscale Scales in Science and Engineering. In: Fish J (ed)
crack modeling. Implementational and computational Multiscale modeling with extended bridging domain
aspects. Int J Numer Methods Eng 89:192–226 method. Oxford University Press, Oxford
51. Nguyen-Xuan H, Hoang T, Nguyen VP (2014) An isoge- 66. Rabczuk T (2013) Computational methods for fracture
ometric analysis for elliptic homogenization problems. in brittle and quasi-brittle solids: state-of-the-art review
Comput Methods Appl Mech Eng 67(9):1722–1741 and future perspectives. Appl Math 2013:849231
52. Svenning E, Fagerström M, Larsson F (2016) On com- 67. Ghayour M, Hosseini-Toudeshky H, Jalalvand M, Bar-
putational homogenization of microscale crack propa- bero EJ (2016) Micro/macro approach for prediction
gation. Int J Numer Meth Eng 108:76–90 of matrix cracking evolution in laminated composites.
53. Zhang R, Zhang L, Wang R, Zhao Y, Huang R (2016) J Compos Mater 50(19):2647–2659
Simulation of a multistage fractured horizontal well 68. Kerfriden P, Passieux JC, Bordas SPA (2012) Local/
with finite conductivity in composite shale gas reser- global model order reduction strategy for the simula-
voir through finite-element method. ACS Energy Fuels tion of quasi-brittle fracture. Int J Numer Methods Eng
30:9036–9049 89:154–179
54. Tene M, Kobaisi MSA, Hajibeygi H (2016) Algebraic 69. Ojo SO, Budarapu PR, Paggi M (2017) A nonlocal adaptive
multiscale method for flow in heterogeneous porous discrete empirical interpolation method combined with
media with embedded discrete fractures (F-AMS). J modified hp-refinement for order reduction of molecular
Comput Phys 321:819–845 dynamics systems. Comput Mater Sci 140:189–208
368
70. Budarapu PR, Reinoso J, Paggi M (2017) Concur- 86. Liang T, Shan T-R, Cheng Y-T, Devine BD, Noordhoek
rently coupled solid shell-based adaptive multiscale M, Li Y, Lu Z, Phillpot SR, Sinnott SB (2013) Classi-
method for fracture. Comput Methods Appl Mech Eng cal atomistic simulations of surfaces and heterogene-
319:338–365. doi:10.1016/j.cma.2017.02.023 ous interfaces with the charge-optimized many body
71. Bischoff M, Ramm E (1997) Shear deformable shell (COMB) potentials. Mater Sci Eng R 74:255–279
elements for large strains and rotations. Int J Numer 87. Ramisetti SB, Anciaux G, Molinari JF (2013) Spatial
Methods Eng 40:4427–4449 filters for bridging molecular dynamics with finite ele-
72. Miehe C (1998) A theoretical and computational model ments at finite temperatures. Comput Methods Appl
for isotropic elastoplastic stress analysis in shells at large Mech Eng 253:28–38
strains. Comput Methods Appl Mech Eng 155:193–233 88. Mulay SS, Becker G, Vayrette R, Raskin JP, Pardoen T,
73. Reinoso J, Paggi M (2014) A consistent interface ele- Galceran M, Godet S, Noels L (2015) Multiscale model-
ment formulation for geometrical and material nonlin- ling framework for the fracture of thin brittle polycrys-
earities. Comput Mech 54(6):1569–1581 talline films: application to polysilicon. Comput Mech
74. Reinoso J, Blázquez A (2016) Application and finite ele- 55:73–91
ment implementation of 7-parameter shell element for 89. Verlet L (1967) Computer “experiments” on classical
geometrically nonlinear analysis of layered CFRP com- fluids. I. Thermodynamical properties of Lennard–
posites. Compos Struct 139:263–273 Jones molecules. Phys Rev 159(5):98–103
75. Reinoso J, Paggi M, Areias PMA (2016) A finite element 90. Abraham FF, Broughton JQ, Bernstein N, Kaxiras E
framework for the interplay between delamination and (1998) Spanning the length scales in dynamic simula-
buckling of rubber-like bi-material systems and stretch- tion. Comput Phys 12:538–546
able electronics. J Eur Ceram Soc 36(9):2371–2382 91. Novoselov KS, Geim AK, Morozov SV, Jiang D, Katsnel-
76. Plimpton S (1995) Fast parallel algorithms for short- son MI, Grigorieva IV, Dubonos SV, Firsov AA (2005)
range molecular dynamics. J Comput Phys 117:1–19 Two-dimensional gas of massless Dirac fermions in
77. Sproul A (2003) Solar cells resources for the secondary graphene. Nature 438:197–200
science teacher, chapter understanding the p-n junc- 92. Khare R, Mielke SL, Paci JT, Zhang S, Ballarini R, Schatz
tion. University of New South Wales, Sydney GC, Belytschko T (2007) Coupled quantum mechani-
78. Köntges M, Kunze I, Kajari-Schröder S, Breitenmoser X, cal/molecular mechanical modeling of the fracture of
Bjrneklett B (2011) The risk of power loss in crystalline defective carbon nanotubes and graphene sheets. Phys
silicon based photovoltaic modules due to microcracks. Rev B 75:075412
Sol Energy Mater Sol Cells 95:1131–1137 93. Lu Q, Gao W, Huang R (2011) Atomistic simulation
79. Paggi M, Sapora A (2013) Numerical modelling of and continuum modeling of graphene nanoribbons
microcracking in PV modules induced by thermo- under uniaxial tension. Modell Simul Mater Sci Eng
mechanical loads. Energy Proc 38:506–515 19:54006
80. Paggi M, Berardone I, Infuso A, Corrado M (2014) 94. Novoselov KS, Falko VI, Colombo L, Gellert PR,
Fatigue degradation and electric recovery in Silicon Schwab MG, Kim K (2012) A roadmap for graphene.
solar cells embedded in photovoltaic modules. Sci Rep Nature 490:192–400
4:4506 95. Ansari R, Ajori S, Motevalli B (2012) Mechanical prop-
81. Käsewieter J, Haase F, Larrodé MH, Köntges M (2014) erties of defective single-layered graphene sheets via
Cracks in solar cell metallization leading to mod- molecular dynamics simulation. Superlattices Micro-
ule power loss under mechanical loads. Energy Proc struct 51:274–289
27:469–477 96. Xu M, Tabarraei A, Paci J, Oswald J, Belytschko T (2012)
82. Yang H, Wang H, Cao D, Sun D, Ju X (2015) Analysis A coupled quantum/continuum mechanics study of
of power loss for crystalline silicon solar module dur- graphene fracture. Int J Fract 173:163–173
ing the course of encapsulation. Int J Photoenergy 97. Budarapu PR, Javvaji B, Sutrakar VK, Mahapatra DR, Zi
2015:251615 G, Rabczuk T (2015) Crack propagation in graphene. J
83. Paggi M, Kajari-Schröder S, Eitner U (2011) Thermo- Appl Phys 118:064307
mechanical deformations in photovoltaic laminates. J 98. Budarapu PR, Javvaji B, Sutrakar VK, Mahapatra DR,
Strain Anal Eng Des 46(8):772–782 Zi G, Paggi M, Rabczuk T (2017) Lattice orientation
84. Infuso A, Corrado M, Paggi M (2014) Image analysis and crack size effect on the mechanical properties
of polycrystalline solar cells and modelling of inter- of graphene. Int J Fract 203(1):81–91. doi:10.1007/
granular and transgranular cracking. J Eur Ceram Soc s10704-016-0115-9
34(11):2713–2722 99. Javvaji B, Budarapu PR, Sutrakar VK, Mahapatra DR, Zi
85. Paggi M, Corrado M, Berardone I (2016) A global/local G, Paggi M, Rabczuk T (2016) Mechanical properties of
approach for the prediction of the electric response of graphene: molecular dynamics simulations correlated
cracked solar cells in photovoltaic modules under the to continuum based scaling laws. Comput Mater Sci
action of mechanical loads. Eng Fract Mech 168:40–57 125:319–327
369
100. Amani J, Oterkus E, Areias P, Zi G, Nguyen-Thoi T, performance molecular simulations through multi-
Rabczuk T (2016) A non-ordinary state-based peridy- level parallelism from laptops to supercomputers. Soft-
namics formulation for thermoplastic fracture. Int J wareX 1–2:19–25
Impact Eng 87:83–94 114. Barkaoui A, Tlili B, Vercher-Martńez A (2016) A mul-
101. Wang HS (2015) A meshfree variational multiscale tiscale modelling of bone ultrastructure elastic pro-
methods for thermo-mechanical material failure. Theo- prieties using finite elements simulation and neural
ret Appl Fract Mech 75:1–7 network method. Comput Methods Programs Biomed
102. Miehe C, Vallicotti D, Zäh DD (2015) Computational 134:69–78
structural and material stability analysis in finite 115. Harvey M, Giupponi G, De Fabritiis G (2009) ACEMD:
electro-elasto-statics of electro-active materials. Int J accelerated molecular dynamics simulations in the
Numer Methods Eng 102:1605–1637 microseconds timescale. J Chem Theory Comput
103. Thomas S, Ajith KM (2014) Molecular dynamics simu- 5:1632
lation of the thermo-mechanical properties of mon- 116. Harvey M, De Fabritiis G (2009) An implementation
olayer graphene sheet. Proc Mater Sci 5:489–498 of the smooth particle-mesh Ewald (PME) method on
104. Chih-Ping W, Kuan-Hao C, Ming-Wang Y (2008) A GPU hardware. J Chem Theory Comput 5:2371–2377
meshfree DRK-based collocation method for the cou- 117. Agilemolecule (2016) http://www.biomolecular-mode-
pled analysis of functionally graded magneto-elec- ling.com/abalone/. Stockholm University
tro-elastic shells and plates. Comput Model Eng Sci 118. Berendsen HJC, van der Spoel D, van Drunen R (1995)
35:181–214 GROMACS: a message-passing parallel molecular
105. Chih-Ping W, Jian-Sin W, Ming-Wang Y (2009) A DRK dynamics implementation. Comput Phys Commun
interpolation-based collocation method for the analysis 91(1):43–56
of functionally graded piezoelectric hollow cylinders 119. Pearlman DA, Case DA, Caldwell JW, Ross WS,
under electro-mechanical loads. Comput Model Eng Cheatham TE, Bolt SD, Ferguson D, Seibel G, Kollman
Sci 52:1–37 P (1995) AMBER, a package of computer programs for
106. Rodríguez GD, Tapia A, Siedel GD, Avilés F (2016) applying molecular mechanics, normal mode analysis,
Influence of structural defects on the electrical proper- molecular dynamics and free energy calculations to
ties of carbon nanotubes and their polymer composites. simulate the structural and energetic properties of mol-
Adv Funct Mater. doi:10.1002/adem.201600116 ecules. Comput Phys Commun 91(1):1–41
107. Tiwary CS, Javvaji B, Kumar C, Mahapatra DR, Ozden 120. Martínez L, Andrade R, Birgin EG, Martínez JM (2009)
S, Ajayan PM, Chattopadhyay K (2015) Chemical-free PACKMOL: a package for building initial configura-
graphene by unzipping carbon nanotubes using cryo- tions for molecular dynamics simulations. J Comput
milling. Carbon. doi:10.1016/j.carbon.2015.03.036 Chem 30(13):2157–2164
108. Tiwary CS, Vishnu D, Kole AK, Brahmanandam J, 121. Humphrey W, Dalke A, Schulten K (1996) VMD: visual
Mahapatra DR, Kumbhakar P, Chattopadhyay K (2015) molecular dynamics. J Mol Graph 14(1):33–38
Stabilization of the high-temperature and high-pres- 122. Qian D, Wagner GJ, Liu WK (2003) A multi-
sure cubic phase of ZnO by temperature-controlled scale projection method for the analysis of car-
milling. J Mater Sci. doi:10.1007/s10853-015-9394-1 bon nanotubes. Comput Methods Appl Mech Eng
109. Vu-Bac N, Lahmer T, Zhang Y, Zhuang X, Rabczuk T 193(17–20):1603–1632
(2014) Stochastic predictions of interfacial character- 123. Qian D, Gondhalekar RH (2004) A virtual atom clus-
istic of polymeric nanocomposites (PNCs). Compos B ter approach to the mechanics of nanostructures. Int J
59:80–95 Multiscale Comput Eng 2(2):277–289
110. Scott WRP, Hünenberger PH, Tironi IG, Mark AE, Bil- 124. Qian D, Chirputkar S (2014) Bridging scale simulation
leter SR, Fennen J, Torda AE, Huber T, Krüger P, van of lattice fracture using enriched space–time finite ele-
Gunsteren WF (1999) The GROMOS biomolecular ment method. Int J Numer Methods Eng 97:819–850
simulation package. J Phys Chem A 103:3596–3607 125. Huang T, Zhang YX, Yang C (2016) Multiscale model-
111. Schmid N, Christ CD, Christen M, Eichenberger AP, ling of multiple-cracking tensile fracture behaviour of
van Gunsteren WF (2012) Architecture, implemen- engineered cementitious composites. Eng Fract Mech
tation and parallelisation of the GROMOS software 160:52–66
for biomolecular simulation. Comput Phys Commun 126. Fereidoon A, Rajabpour M, Hemmatian H (2013)
183:890–903 Fracture analysis of epoxy/SWCNT nanocomposite
112. Ferreira RJ, Ferreira MJU, dos Santos DJVA (2012) based on global–local finite element model. Compos B
Insights on p-glycoprotein’s efflux mechanism obtained 54:400–408
by molecular dynamics simulations. J Chem Theory 127. Melenk JM, Babuska I (1996) The partition of unity
Comput 8(6):1853–1864 finite element method: basic theory and applications.
113. Abraham MJ, Murtola T, Schulz R, Páll S, Smith Comput Methods Appl Mech Eng 139(1–4):289–314
JC, Hess B, Lindahl E (2015) GROMACS: high
370
128. Cai Y, Zhuang X, Augarde C (2010) A new partition of mechanics problems with singular field of arbitrary
unity finite element free from linear dependence prob- order. Comput Methods Appl Mech Eng 253:252–273
lem and processing delta property. Comput Methods 144. Vu-Bac N, Nguyen-Xuan H, Chen L, Bordas SPA, Ker-
Appl Mech Eng 199:1036–1043 friden P, Simpson RN, Liu GR, Rabczuk T (2011) A
129. Rabczuk T, Belytschko T, Xiao SP (2004) Stable particle noded-based smoothed XFEM for fracture mechanics.
methods based on Lagrangian kernels. Comput Meth- Comput Model Eng Sci 73:331–356
ods Appl Mech Eng 193(12–14):1035–1063 145. Strouboulis T, Copps K, Babuska I (2000) The gener-
130. Belytschko T, Lu YY, Gu L (1994) Element-free galerkin alized finite element method: an example of its imple-
methods. Int J Numer Methods Eng 37:229–256 mentation and illustration of its performance. Int J
131. Liu WK, Jun S, Zhang YF (1995) Reproducing ker- Numer Methods Eng 47:1401–1417
nel particle methods. Int J Numer Methods Fluids 146. Strouboulis T, Copps K, Babuska I (2001) The gener-
20:1081–1106 alized finite element method. Comput Methods Appl
132. Liu WK, Jun S, Li S, Adee J, Belytschko T (1995b) Mech Eng 190(32–33):4081–4193
Reproducing kernel particle methods for structural 147. Strouboulis T, Zhang L, Babuska I (2003) Generalized
dynamics. Int J Numer Methods Eng 38:1655–1679 finite element method using mesh-based handbooks:
133. Daux CC, Möes N, Dolbow J, Sukumar N, Belytschko T application to problems in domains with many voids.
(2000) Arbitrary branched and intersecting cracks with Comput Methods Appl Mech Eng 192(28):3109–3161
the extended finite element method. Int J Numer Meth- 148. Strouboulis T, Babuska I, Hidajat R (2006) The gener-
ods Eng 48:1741–1760 alized finite element method for Helmholtz equation:
134. Sukumar N, Moës N, Moran B, Belytschko T (2000) theory, computation, and open problems. Comput
Extended finite element method for three-dimen- Methods Appl Mech Eng 195:4711–4731
sional crack modeling. Int J Numer Methods Eng 149. Strouboulis T, Hidajat R, Babuska I (2008) The general-
48:1549–1570 ized finite element method for helmholtz equation, Part
135. Zi G, Chen H, Xu JX, Belytschko T (2005) The extended II: effect of choice of handbook functions, error due
finite element method for dynamic fractures. Shock Vib to absorbing boundary conditions and its assessment.
12(1):9–23 Comput Methods Appl Mech Eng 197:364–380
136. Areias PMA, Belytschko T (2005) Non-linear analysis 150. Babuska I, Banerjee U, Osborn JE (2004) Generalized
of shells with arbitrary evolving cracks using XFEM. Int finite element methods: main ideas, results, and per-
J Numer Methods Eng 62:384–415 spective. Int J Comput Methods 1(1):67–103
137. Areias PMA, Belytschko T (2005) Analysis of three- 151. Rabczuk T, Zi G (2007a) A meshfree method based on
dimensional crack initiation and propagation using the the local partition of unity for cohesive cracks. Comput
extended finite element method. Int J Numer Methods Mech 39:743–760
Eng 63:760–788 152. Rabczuk T, Bordas SP, Zi G (2007b) A three-dimen-
138. Bordas SPA, Rabczuk T, Nguyen-Xuan H, Natara- sional meshfree method for continuous multiple-crack
jan S, Bog T, Nguyen VP, Do MQ, Nguyen VH (2010) initiation, propagation and junction in statics and
Strain smoothing in FEM and XFEM. Comput Struct dynamics. Comput Mech 40(3):473–495
88(23–24):1419–1443 153. Zi G, Rabczuk T, Wall W (2007) Extended meshfree
139. Bordas SPA, Natarajan S, Pont SD, Rabczuk T, Kerfriden methods without branch enrichment for cohesive
P, Mahapatra DR, Noel D, Gao Z, Gao Z (2011) On the cracks. Comput Mech 40(2):367–382
performance of strain smoothing for enriched finite 154. Bordas SPA, Rabczuk T, Zi G (2008) Three-dimensional
element approximations (XFEM/GFEM/PUFEM). Int J crack initiation, propagation, branching and junc-
Numer Methods Eng 86(4–5):637–666 tion in non-linear materials by an extended meshfree
140. Nanthakumar SS, Lahmer T, Rabczuk T (2013) Detec- method without asymptotic enrichment. Eng Fract
tion of flaws in piezoelectric structures using xfem. Int J Mech 75(5):943–960
Numer Methods Eng 96(6):373–389 155. Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H (2008) A
141. Nanthakumar SS, Lahmer T, Rabczuk T (2014) Detec- geometrically non-linear three dimensional cohesive
tion of multiple flaws in piezoelectric structures using crack method for reinforced concrete structures. Eng
XFEM and level sets. Comput Methods Appl Mech Eng Fract Mech 75(16):4740–4758
275:98112 156. Rabczuk T, Samaniego E (2008) Discontinuous model-
142. Nguyen-Thanh N, Kiendl J, Nguyen-Xuan H, Wuch- ling of shear bands using adaptive meshfree methods.
ner R, Bletzinger KU, Bazilevs Y, Rabczuk T (2011) Comput Methods Appl Mech Eng 197(6–8):641–658
Rotation free isogeometric thin shell analysis using 157. Rabczuk T, Gracie R, Song J-H, Belytschko T (2010)
PHT-splines. Comput Methods Appl Mech Eng Immersed particle method for fluid-structure interac-
200(47–48):3410–3424 tion. Int J Numer Methods Eng 81(1):48–71
143. Nguyen-Xuan H, Liu GR, Bordas SPA, Natarajan S, 158. Ventura G, Xu XJ, Belytschko T (2002) A vector level
Rabczuk T (2013) An adaptive singular ES-FEM for set method and new discontinuity approximations
371
for crack growth by EFG. Int J Numer Methods Eng 174. Guowei M, Xinmei A, Lei H (2010) The numerical
54:923–944 manifold method: a review. Int J Comput Methods
159. Rabczuk T, Belytschko T (2004) Cracking particles: 7(1):1–32
a simplified meshfree method for arbitrary evolving 175. Miehe C, Welschinger F, Hofacker M (2010) Thermo-
cracks. Int J Numer Methods Eng 61:2316–2343 dynamically consistent phase-field models of fracture:
160. Rabczuk T, Areias PMA (2006) A new approach variational principles and multi-field FE implementa-
for modelling slip lines in geological materials with tions. Int J Numer Methods Eng 83(10):1273–1311
cohesive models. Int J Numer Anal Methods Eng 176. Miehe C, Hofacker M, Welschinger F (2010) A phase
30(11):1159–1172 field model for rate-independent crack propagation:
161. Rabczuk T, Belytschko T (2006) Application of particle robust algorithmic implementation based on operator
methods to static fracture of reinforced concrete struc- splits. Comput Methods Appl Mech Eng 199:2765–2778
tures. Int J Fract 137(1–4):19–49 177. Msekh MA, Sargado M, Jamshidian M, Areias P, Rab-
162. Rabczuk T, Belytschko T (2007) A three dimensional czuk T (2015) ABAQUS implementation of phase-
large deformation meshfree method for arbitrary field model for brittle fracture. Comput Mater Sci
evolving cracks. Comput Methods Appl Mech Eng 96B:472–484
1960(29–30):2777–2799 178. Amiri F, Millán D, Shen Y, Rabczuk T, Arroyo M (2014)
163. Rabczuk T, Areias PMA, Belytschko T (2007) A mesh- Phase-field modeling of fracture in linear thin shells.
free thin shell method for nonlinear dynamic fracture. Theoret Appl Fract Mech 69:102–109
Int J Numer Methods Eng 72(5):524–548 179. Seleson P, Beneddine S, Prudhomme S (2013) A force-
164. Rabczuk T, Song JH, Belytschko T (2009) Simulations based coupling scheme for peridynamics and classical
of instability in dynamic fracture by the cracking parti- elasticity. Comput Mater Sci 66:34–49
cles method. Eng Fract Mech 76(6):730–741 180. Tong Q, Li S (2016) Multiscale coupling of molecu-
165. Rabczuk T, Zi G, Bordas SP, Nguyen-Xuan H (2010) A lar dynamics and peridynamics. J Mech Phys Solids
simple and robust three-dimensional cracking-particle 95:169–187
method without enrichment. Comput Methods Appl 181. Hughes TJR, Cottrell J, Bazilevs Y (2005) Isogeometric
Mech Eng 199(37–40):2437–2455 analysis: CAD, finite elements, NURBS, exact geometry
166. Vu-Bac N, Nguyen-Xuan H, Chen L, Bordas SPA, and mesh refinement. Comput Methods Appl Mech
Zhuang X, Liu GR, Rabczuk T (2013) A phantom- Eng 194:4135–4195
node method with edge-based strain smoothing 182. Bazilevs Y, Calo VM, Cottrell JA, Evans JA, Hughes TJR,
for linear elastic fracture mechanics. J Appl Math. Lipton S, Scott MA, Sederberg TW (2010) Isogeometric
doi:10.1155/2013/978026 analysis using T-splines. Comput Methods Appl Mech
167. Chau-Dinh T, Zi G, Lee PS, Song JH, Rabczuk T (2010) Eng 199(5–8):229–263
Phantom-node method for shell models with arbitrary 183. Nguyen VP, Nguyen-Xuan H (2013) High-order
cracks. Comput Struct 92–93:242–256 b-splines based finite elements for the delamina-
168. Rabczuk T, Zi G, Gerstenberger A, Wall WA (2008) A tion analysis of laminated composites. Compos Struct
new crack tip element for the phantom node method 102:261–275
with arbitrary cohesive cracks. Int J Numer Methods 184. Chien HT, Ferreira AJ, Carrera E, Nguyen-Xuan H
Eng 75(5):577–599 (2013) Isogeometric analysis of laminated composite
169. Song JH, Areias PMA, Belytschko T (2006) A method and sandwich plates using a layerwise deformation the-
for dynamic crack and shear band propagation ory. Compos Struct 104:196–214
with phantom nodes. Int J Numer Methods Eng 185. Chien HT, Ferreira AJM, Rabczuk T, Bordas SPA,
67(6):868–893 Nguyen-Xuan H (2014) Isogeometric analysis of lami-
170. Hansbo A, Hansbo P (2004) A finite element method nated composite and sandwich plates using a new
for the simulation of strong and weak discontinuities inverse trigonometric shear deformation theory. Eur J
in solid mechanics. Comput Methods Appl Mech Eng Mech A Solids 43:89–108
193(33–35):3523–3540 186. Nguyen-Xuan H, Thai CH, Nguyen-Thoi T (2013)
171. Mergheim J, Kuhl E, Steinmann P (2005) A finite ele- Isogeometric finite element analysis of composite sand-
ment method for the computational modelling of cohe- wich plates using a new higher order shear deformation
sive cracks. Int J Numer Methods Eng 63:276–289 theory. Compos Part B 55:558–574
172. Mergheim J, Kuhl E, Steinmann P (2007) Towards the 187. Areias PMA, Rabczuk T (2013) Finite strain fracture of
algorithmic treatment of 3d strong discontinuities. plates and shells with configurational forces and edge
Commun Numer Methods Eng 23:97–108 rotations. Int J Numer Methods Eng 94:1099–1122
173. Cai Y, Zhuang X, Zhu H (2013) A generalized and 188. Rabczuk T, Areias PMA (2006) A meshfree thin shell
efficient method for finite cover generation in the for arbitrary evolving cracks based on an extrinsic basis.
numerical manifold method. Int J Comput Methods Comput Model Eng Sci 16(2):115–130
10(5):1350028
372
189. Nguyen-Thanh N, Valizadeh N, Nguyen MN, Nguyen- 205. Liu GR, Gu YT (2001) A local radial point interpolation
Xuan H, Zhuang X, Areias P, Zi G, Bazilevs Y, De Loren- method (LRPIM) for free vibration analyses of 2-D sol-
zis L, Rabczuk T (2015) An extended isogeometric thin ids. J Sound Vib 246(1):29–46
shell analysis based on Kirchhoff–Love theory. Comput 206. Liu GR, Zhang GY, Gu YT, Wang YY (2005) A meshfree
Methods Appl Mech Eng 284:265–291 radial point interpolation method (RPIM) for three-
190. Plews JA, Duarte CA (2016) A two-scale generalized dimensional solids. Comput Mech 36(6):421–430
finite element approach for modeling localized thermo- 207. Shih-Wei Y, Yung-Ming W, Chih-Ping W, Hsuan-Teh
plasticity. Int J Numer Methods Eng 108:1123–1158 H (2010) A meshless collocation method based on the
191. Akkutlu IY, Efendiev Y, Vasilyeva M (2016) Multiscale differential reproducing kernel approximation. Comput
model reduction for shale gas transport in fractured Model Eng Sci 60:1–39
media. Comput Geosci 20:953–973 208. Zhuang X, Augarde C, Mathisen K (2012) Fracture
192. Efendiev E (2015) Generalized multiscale finite element modelling using meshless methods and level sets in 3D:
methods (GMsFEM). J Comput Phys 251:116–135 framework and modelling. Int J Numer Methods Eng
193. Belytschko T, Black T (1999) Elastic crack growth in 92:969–998
finite elements with minimal remeshing. Int J Numer 209. Chih-Ping W, Shih-Wei Y, Ming-Wang Y, Hsuan-Teh
Methods Eng 45(5):601–620 H (2011) A meshless collocation method for the plane
194. Moës N, Dolbow J, Belytschko T (1999) A finite ele- problems of functionally graded material beams and
ment method for crack growth without remeshing. Int J plates using the DRK interpolation. Mech Res Com-
Numer Methods Eng 46(1):131–150 mun 38:471–476
195. Belytschko T, Gracie R, Ventura G (2009) A review of 210. Yung-Ming W, Syuan-Mu C, Chih-Ping W (2010) A
extended/generalized finite element methods for mate- meshless collocation method based on the differen-
rial modeling. Model Simul Mater Sci Eng 17:043001 tial reproducing kernel interpolation. Comput Mech
196. Fries T-P, Belytschko T (2010) The extended/gen- 45:585–606
eralized finite element method: an overview of the 211. Griffith AA (1921) The phenomena of rapture and flow
method and its applications. Int J Numer Methods Eng in solids. Philos Trans R Soc Lond A 221:163–198
84(3):253–304 212. Irwin G (1957) Analysis of stresses and strains near
197. Karihaloo BL, Xiao QZ (2003) Modelling of stationary the end of a crack traversing a plate. J Appl Mech
and growing cracks in fe framework without remeshing: 24:361–364
a state-of the-art review. Comput Struct 81(3):119–129 213. Kuhn C, Müller R (2010) A continuum phase field
198. Mohammadi S (2008) Extended finite element method model for fracture. Eng Fract Mech 77(18):3625–3634
for fracture analysis of structures, vol Oxford. Blackwell 214. Karma A, Kessler DA, Levine H (2001) Phase-field
Publishing, Hoboken model of mode III dynamic fracture. Phys Rev Lett
199. Nguyen VP, Rabczuk T, Bordas SP, Duflot M (2008) 87:045501
Meshless methods: a review and computer implementa- 215. Borden MJ, Hughes TJR, Landis CM, Verhoosel CV
tion aspects. Math Comput Simul 79(3):763–813 (2014) A higher-order phase-field model for brit-
200. Belytschko T, Tabbara M (1996) Dynamic fracture tle fracture. Comput Methods Appl Mech Eng
using element-free Galerkin methods. Int J Numer 273:100–118
Methods Eng 39:923–938 216. Areias P, Rabczuk T, Msekh MA (2016) Phase-field
201. Bauman P, Dhia H, Elkhodja N, Oden J, Prudhomme analysis of finite-strain plates and shells including ele-
S (2008) On the application of the Arlequin method to ment subdivision. Comput Methods Appl Mech Eng
the coupling of particle and continuum models. Com- 312:322–350
put Mech 42(4):511–530 217. Yingjun G, Zhirong L, Lilin H, Hong M (2016) Phase
202. Amiri F, Anitescu C, Arroyo M, Bordas S, Rabczuk T field crystal study of nano-crack growth and branch in
(2014) XLME interpolants, a seamless bridge between materials. Modell Simul Mater Sci Eng 24:055010
xfem and enriched meshless methods. Comput Mech 218. Giovanardi B, Scotti A, Formaggia L (2017) A hybrid
53(1):45–57 XFEM-phase field (Xfield) method for crack propaga-
203. Atluri SN, Shengping S (2002) The meshless local tion in brittle elastic materials. Comput Methods Appl
Petrov-Galerkin (MLPG) method: a simple and less- Mech Eng. doi:10.1016/j.cma.2017.03.039
costly alternative to the finite element and boundary 219. Areias PMA, Song JH, Belytschko T (2006) Analysis of
element methods. Comput Model Eng Sci 3(1):11–51 fracture in thin shells by overlapping paired elements.
204. Liu K, Long S, Li G (2006) A simple and less-costly Comput Methods Appl Mech Eng 195:5343–5360
meshless local Petrov-Galerkin (MLPG) method for the 220. Liu WK, Qian D, Gonella S, Li S, Chen W, Chirputkar
dynamic fracture problem. Eng Anal Boundary Elem S (2010) Multiscale methods for mechanical science of
30(1):72–76 complex materials: bridging from quantum to stochas-
tic multiresolution continuum. Int J Numer Methods
Eng 83:1039–1080
373
221. Chen H, Zang M, Zhang YX (2016) A ghost particle- 236. Gurtin ME, Podio-Guidugli P (1996) Configurational
based coupling approach for the combined finite-dis- forces and the basic laws for crack propagation. J Mech
crete element method. Finite Elem Anal Des 114:68–77 Phys Solids 44(6):905–927
222. Kojic M, Filipovic N, Tsuda A (2008) A mesoscopic 237. Gurtin ME, Podio-Guidugli P (1998) Configurational
bridging scale method for fluids and coupling dis- forces and a constitutive theory for crack propagation
sipative particle dynamics with continuum finite ele- that allows for kinking and curving. J Mech Phys Solids
ment method. Comput Methods Appl Mech Eng 46(8):1343–1378
197(6–8):821–833 238. Miehe C, Gürses E (2007) A robust algorithm for con-
223. Li X, Wan K (2011) A bridging scale method for granu- figurational-force-driven brittle crack propagation with
lar materials with discrete particle assembly-cosserat radaptive mesh alignment. Int J Numer Methods Eng
continuum modeling. Comput Geotech 38:1052–1068 72(2):127–157
224. Kadowaki H, Liu WK (2004) Bridging multi-scale 239. Mosler J (2009) A variationally consistent approach
method for localization problems. Comput Methods for crack propagation based on configurational forces.
Appl Mech Eng 193(30–32):3267–3302 In: IUTAM symposium on progress in the theory and
225. Tang S, Hou TY, Liu WK (2006) A mathematical frame- numerics of configurational mechanics, vol 17. IUTAM
work of the bridging scale method. Int J Numer Meth- Bookseries, pp 239–247
ods Eng 65:1688–1713 240. Sih GC (1974) Strain-energy-density factor applied to
226. Wang W, Chang K-H (2013) Continuum-based sen- mixed mode crack problems. Int J Fract 10(3):305–321
sitivity analysis for coupled atomistic and continuum 241. Wu CH (1978) Fracture under combined loads by max-
simulations for 2-d applications using bridging scale imum energy release rate criterion. J Appl Mech Trans
decomposition. Struct Multidiscipl Optim 47:867–892 ASME 45(3):553–558
227. Guidault PA, Belytschko T (2007) On the L2 and the H 1 242. Goldstein RV, Salganik RL (1974) Brittle fracture of sol-
couplings for an overlapping domain decomposition ids with arbitrary cracks. Int J Fract 10(4):507–523
method using lagrange multipliers. Int J Numer Meth- 243. Benedetti I, Aliabadi MH (2015) Multiscale modeling of
ods Eng 70:322–350 polycrystalline materials: a boundary element approach
228. Gracie R, Belytschko T (2011) Adaptive continuum- to material degradation and fracture. Comput Methods
atomistic simulations of dislocation dynamics. Int J Appl Mech Eng 289:429–453
Numer Methods Eng 86(4–5):575–597 244. Cynthia LK, Plimpton SJ, Hamilton JC (1998) Dislo-
229. Belytschko T, Xiao SP (2004) A bridging domain cation nucleation and defect structure during surface
method for coupling continua with molecu- indentation. Phys Rev B 58:11085
lar dynamics. Comput Methods Appl Mech Eng 245. Budarapu PR, Gracie R, Yang S-W, Zhuang X, Rabczuk
193(17–20):1645–1669 T (2014) Efficient coarse graining in multiscale mod-
230. Fish J, Chen W (2004) Discrete-to-continuum bridging eling of fracture. Theoret Appl Fract Mech 69:126–143
based on multigrid principles. Comput Methods Appl 246. Dipasquale D, Zaccariotto M, Galvanetto U (2014)
Mech Eng 193:1693–1711 Crack propagation with adaptive grid refinement in 2D
231. Tabarraei A, Wang X, Sadeghirad A, Song JH (2014) An peridynamics. Int J Fract 190:1–22
enhanced bridging domain method for linking atom- 247. Unger JF (2013) An fe 2-x 1 approach for multi-
istic and continuum domains. Finite Elem Anal Des scale localization phenomena. J Mech Phys Solids
92:36–49 61:928–948
232. Xu M, Belytschko T (2008) Conservation properties 248. Park K, Paulino GH, Celes W, Espinha R (2012) Adap-
of the bridging domain method for coupled molecu- tive mesh refinement and coarsening for cohesive zone
lar/continuum dynamics. Int J Numer Methods Eng modeling of dynamic fracture. Int J Numer Methods
76:278–294 Eng 92:1–35
233. Xu M, Gracie R, Belytschko T (2010) A continuum-to- 249. Belytschko T, Song JH (2010) Coarse-graining of mul-
atomistic bridging domain method for composite lat- tiscale crack propagation. Int J Numer Methods Eng
tices. Int J Numer Methods Eng 81:1635–1658 81(5):537–563
234. Anciaux G, Ramisetti SB, Molinari JF (2012) A finite 250. Miller RE, Tadmor EB (2009) A unified framework and
temperature bridging domain method for MD-FE cou- performance benchmark of fourteen multiscale atom-
pling and application to a contact problem. Comput istic/continuum coupling methods. Model Simul Mater
Methods Appl Mech Eng 205–208:204–212 Sci Eng 17:053001
235. Tu F, Ling D, Bu L, Yang Q (2014) Generalized bridg- 251. Nair AK, Warner DH, Hennig RG (2011) Coupled
ing domain method for coupling finite elements with quantum-continuum analysis of crack tip processes in
discrete elements. Comput Methods Appl Mech Eng aluminum. J Mech Phys Solids 59(12):2476–2487
276:509–533 252. Liu W, Hong J-W (2012) A coupling approach of dis-
cretized peridynamics with finite element method.
Comput Methods Appl Mech Eng 245–246:163–175
374
253. Talebian M, Al-Khoury R, Sluys LJ (2013) A computa- 266. Silling SA (2000) Reformulation of elasticity theory for
tional model for coupled multiphysics processes of CO2 discontinuities and long-range forces. J Mech Phys Sol-
sequestration in fractured porous media. Adv Water ids 48:175–209
Resour 59:238–255 267. Silling SA, Lehoucq RB (2010) Peridynamic theory of
254. Jothi S, Croft TN, Brown SGR (2015) Multiscale mul- solid mechanics. Adv Appl Mech 44:73–166
tiphysics model for hydrogen embrittlement in poly- 268. Ren H, Zhuang X, Cai Y, Rabczuk T (2016) Dual-
crystalline nickel. J Alloys Compd 645:S500–S504 horizon peridynamics. Int J Numer Meth Eng
255. Reifsnider K, Raihan MDR, Vadlamudi V (2016) Heter- 108(12):1451–1476
ogeneous fracture mechanics for multi-defect analysis. 269. Ren H, Zhuang X, Rabczuk T (2017) Dual-horizon
Compos Struct 156:20–28 peridynamics: A stable solution to varying horizons.
256. Grujicic M, Ramaswami S, Yavari R, Galgalikar R Comput Methods Appl Mech Eng 318:762–782
(2016) Multiphysics computational analysis of white- 270. Parks ML, Lehoucq RB, Plimpton SJ, Silling SA (2008)
etch cracking failure mode in wind turbine gear- Implementing peridynamics within a molecular
box bearings. Proc IMechE Part L J Mater Des Appl dynamics code. Comput Phys Commun 179:777–783
230(1):43–63 271. Rahman R, Foster JT (2014) Bridging the length scales
257. Wu W, Xiao X, Huang X, Yan S (2014) A multiphysics through nonlocal hierarchical multiscale modeling
model for the in situ stress analysis of the separator in a scheme. Comput Mater Sci 92:401–415
lithium-ion battery cell. Comput Mater Sci 83:127–136 272. Brothers MD, Foster JT, Millwater HR (2014) A com-
258. Jiang T, Rudraraju S, Roy A, Van der Ven A, Garikipati parison of different methods for calculating tangent-
K, Falk ML (2016) Multiphysics simulations of lithia- stiffness matrices in a massively parallel computational
tion-induced stress in Li1+x Ti2O 4 electrode particles. J peridynamics code. Comput Methods Appl Mech Eng
Phys Chem C 120:27871–27881 279:247–267
259. PantanoMF Bernai RA, Pagnotta L, Espinosa HD 273. Zhang Y, Yao H, Ortiz C, Xu J, Dao M (2012) Bio-
(2015) Multiphysics design and implementation of a inspired interfacial strengthening strategy through geo-
microsystem for displacement-controlled tensile testing metrically interlocking designs. J Mech Behav Biomed
of nanomaterials. Meccanica 50:549–569 Mater 15:70–77
260. Saft A, Kaliske M (2013) A hybrid interface-element for 274. Egan P, Sinko R, LeDuc PR, Keten S (2005) The role of
the simulation of moisture-induced cracks in wood. mechanics in biological and bio-inspired systems. Nat
Eng Fract Mech 102:32–50 Commun 6:7418
261. Yu M-F, Lourie O, Dyer MJ, Moloni K, Kelly TF, Ruoff 275. Tang Z, Kotov NA, Magonov S, Ozturk B (2003) Nano-
RS (2000) Strength and breaking mechanism of mul- structured artificial nacre. Nat Mater 2(6):413
tiwalled carbon nanotubes under tensile load. Science 276. Niebel TP, Bouville F, Kokkinis D, Studart AR (2016)
287:637–640 Role of the polymer phase in the mechanics of nacre-
262. Mielke SL, Troya D, Zhang S, Li J-L, Xiao S, Car R, Ruoff like composites. J Mech Phys Solids 96:133–146
RS, Schatz GC, Belytschko T (2004) The role of vacancy 277. Awaja F, Zhang S, Tripathi M, Nikiforov A, Pugno N
defects and holes in the fracture of carbon nanotubes. (2016) Cracks, microcracks and fracture in polymer
Chem Phys Lett 390:413–420 structures: formation, detection, autonomic repair.
263. Svensson MM, Humbel S, Froese RDJ, Matsubara T, Prog Mater Sci 83:536–573
Sieber S, Morokuma K (1996) ONIOM: a multilayered 278. Wang W, Elbanna A (2014) Crack propagation in bone
integrated MO + MM method for geometry optimi- on the scale of mineralized collagen fibrils: role of poly-
zations and single point energy predictions. A test for mers with sacrificial bonds and hidden length. Bone
Diels Alder reactions and Pt(P(t-Bu)3)2+H2 oxidative 68:20–31
addition. J Phys Chem 100:19357–19363 279. Fratzl P, Kolednik O, Fischer FD, Dean MN (2016) The
264. Khare R, Mielke SL, Schatz GC, Belytschko T (2008) mechanics of tessellations-bioinspired strategies for
Multiscale coupling schemes spanning the quan- fracture resistance. Chem Soc Rev 45:252
tum mechanical, atomistic forcefield, and con- 280. Ma S, Scheider I, Bargmann S (2016) Continuum dam-
tinuum regimes. Comput Methods Appl Mech Eng age modeling and simulation of hierarchical dental
197:3190–3202 enamel. Model Simul Mater Sci Eng 24:045014
265. Park JY, Park CH, Park JS, Kong K-J, Chang H, Im S 281. Ural A, Mischinski S (2013) Multiscale modeling of
(2010) Multiscale computations for carbon nanotubes bone fracture using cohesive finite elements. Eng Fract
based on a hybrid QM/QC (quantum mechanical Mech 103:141–152
and quasicontinuum) approach. J Mech Phys Solids 282. Rabiei R, Bekah S, Barthelat F (2010) Failure mode
58:86–102 transition in nacre and bone-like materials. Acta Bio-
mater 6:4081–4089
375
Dr. Budarapu is a faculty member in the is Computational Solid Mechanics with emphasis on
School of Mechanical Sciences, Indian method development for problems involving fracture and
Institute of Technology, Bhubaneswar, failure of solids and fluid-structure interaction. He is par-
India. Dr. Budarapu’s research focus is ticularly interested in developing multiscale methods and
mainly on computational methods for frac- their application to computational materials design. His
ture and failure of solids. His primary main research topics include: Constitutive Modeling, Mate-
interests are in developing adaptive multiscale methods and rial Instabilities, Fracture, Strain Localization, Numerical
their application to materials and design for various appli- Methods (Extended Finite Element and Mesfhree Meth-
cations. His areas of interests include: multiscale methods ods), Isogeometric Analysis, Computational Fluid-Struc-
for fracture, molecular dynamics, fracture studies in mul- ture Interaction and Biomechanical Engineering. He is the
tiphysics problems and polymer nano-composites. ISI Highly Cited Researcher (2014–2016) in the categories
‘Computer Science’ and ‘Engineering’.
T. Rabczuk is professor of computational
mechanics at Institut für Strukturmechanik
(ISM), Bauhaus Universität-Weimar, since
2009. The research focus of Prof. Rabczuk
376

20-Multiscale Methods For A Fracture A Review

Uploaded by

Copyright:

Available Formats

20-Multiscale Methods For A Fracture A Review

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

20-Multiscale Methods For A Fracture A Review

Uploaded by

Copyright:

Available Formats

J. Indian Inst. Sci.

A Multidisciplinary Reviews Journal

© Indian Institute of Science 2017.

Abstract | The global response of a system is often governed by the

1 Introduction to reveal the fundamental mechanics of material

J. Indian Inst. Sci. | VOL 97:3 | 339–376 September 2017| journal.iisc.ernet.in 13

(a) Nano (c)

coarse fine coarse fine coarse fine

Figure 1: Schematics of a hierarchical, b semi-concurrent, and c concurrent multiscale methods. In the

E tot (q, r) = E es (q, r) + E short (q, r)

(a) (b) (c)

3.4 Phase‑field method for fracture

 for further details.

integrated along the crack path Ŵc, which can be in which

(a) straight crack p 0

(b) angled crack

phantom node real node

(b) Solid shell element

crack surface solid shell 1 2

(phantom nodes) coupling (BSM)

Ghost atom Fine scale

Figure 7: a Schematic of a three-dimensional coupled continuum-atomistic model. b Mechanics of

the consistent constraint form conserves linear

and a cracking part based on additional degrees

• Calculate the regions to be refined, Pn+1

Area to be refined Crack

Area to be coarsened Crack

illustrated in the schematic Fig. 13. As a rule of 6 Computer implementation

(a) (b) (d)

Figure 12: Schematic of meshless equivalent coarse-scale model. a Meshfree particles superimposed

6.1.1 Implementation in LAMMPS the current step is represented by ‘ubary’(=0.05)

Figure 17: Crack propagation of a double-edge notched specimen. a Deformed configuration after 28

materials259, bio applications260, to name a few. 7.1 Quantum mechanics and molecular

extremely challenging261. Hence, as an alternative approach by concurrently combining the QM,

No. Address Title Comments

experimental and numerical studies to under- interdisciplinary collaboration among differ-

Mechanisms of energy dissipation in struc-

You might also like

3.4 Phase‑field method for fracture

for further details.

illustrated in the schematic Fig. 13. As a rule of 6 Computer implementation

6.1.1 Implementation in LAMMPS the current step is represented by ‘ubary’(=0.05)

materials259, bio applications260, to name a few. 7.1 Quantum mechanics and molecular