# Dual-horizon Peridynamics

Huilong Ren  Xiaoying Zhuang<EMAIL_ADDRESS>Yongchang Cai
Timon Rabczuk<EMAIL_ADDRESS>State Key Laboratory of Disaster
Reduction in Civil Engineering, College of Civil Engineering,Tongji
University, Shanghai 200092, China  Institute of Structural Mechanics,
Bauhaus-Universität Weimar, Germany  School of Civil, Environmental and
Architectural Engineering, Korea University, South Korea

###### Abstract

In this paper we develop a new Peridynamic approach that naturally includes
varying horizon sizes and completely solves the ”ghost force” issue.
Therefore, the concept of dual-horizon is introduced to consider the
unbalanced interactions between the particles with different horizon sizes.
The present formulation is proved to fulfill both the balances of linear
momentum and angular momentum. Neither the ”partial stress tensor” nor the
”‘slice” technique are needed to ameliorate the ghost force issue in [1]. The
consistency of reaction forces is naturally fulfilled by a unified simple
formulation. The method can be easily implemented to any existing peridynamics
code with minimal changes. A simple adaptive refinement procedure is proposed
minimizing the computational cost. The method is applied here to the three
Peridynamic formulations, namely bond based, ordinary state based and non-
ordinary state based Peridynamics. Both two- and three- dimensional examples
including the Kalthof-Winkler experiment and plate with branching cracks are
tested to demonstrate the capability of the method in solving wave
propagation, fracture and adaptive analysis .

###### keywords:

Peridynamics, horizon variable, dual-horizon, ghost force, spurious wave
reflection, adaptive refinement

##  1 Introduction

Peridynamics has recently attracted wide interests for researchers in
computational solid mechanics since it provides the possibility to model
dynamic fracture with ease. The crack is part of the solution from PD
simulation instead of part of the problem and no representation of the crack
topology is needed. The original PD method was proposed by Silling [2] in 2000
and has been exploited onwards for extensive applications of mechanical
problems including impact loading, fragmentation [3, 4], composites
delamination [5], beam and plate structures [6, 7].

In PD, the classical balance equations are formulated in an integral form
instead of partial differential form. The particles interact with each other
when the distances between the particles are within a threshold value, called
horizon. The equations of motion at any time $t$ is expressed as

$\displaystyle\rho\mathbf{\ddot{u}}(\mathbf{x},t)=\int_{H_{\mathbf{x}}}\mathbf{f}\left(\mathbf{u}(\mathbf{x}^{\prime},t)-\mathbf{u}(\mathbf{x},t),\mathbf{x}^{\prime}-\mathbf{x}\right)-\mathbf{f}\left(\mathbf{u}(\mathbf{x},t)-\mathbf{u}(\mathbf{x}^{\prime},t),\mathbf{x}-\mathbf{x}^{\prime}\right)\mathbf{d}V_{\mathbf{x}^{\prime}}+\mathbf{b}(\mathbf{x},t)\
,$ (1)

where $H_{\mathbf{x}}$ denotes the horizon (spherical domain) belonging to
$\mathbf{x}$, $\mathbf{u}$ is the displacement vector, $\mathbf{b}$ denotes
the body force, $\rho$ is mass density in the reference configuration, and
$\mathbf{f}$ is a pairwise force function that computes the force vector (per
unit volume squared) [8]. Eq.(1) shows the key idea of PD to unify continuous
and discontinuous media within a single consistent set of equations. The
governing equation is written in an integral form instead of an partial
differential form. The crack surfaces are formed as the outcome of motion and
constitutive models, and there is no entanglement of additional crack
kinematics or geometry treatment. The fracture behaviour including crack
branching and coalescence of multiple cracks is captured through the breakage
of the bonds between particles. Therefore, the need of smoothing of crack
surfaces or branching criterion in the extended finite element method (XFEM)
[9], meshless methods [10] or other partition of unity methods (PUM) [11] is
completely removed. The extension of PD from 2D to 3D problems is greatly
facilitated for the computer implementation, which is not always the case in
other methods [12, 13].

The original PD started with the bond-based formulation (BB-PD) where the
bonds behave like springs and independent of each other. BB-PD can be regarded
as a special case of a more general theory, the state-based peridynamics (SB-
PD) [2, 14] which can be suited for, theoretically, any type of constitutive
model and large deformation analysis. The key difference between SB-PD and BB-
PD is that in the former, the bond deformation depends on collective
deformation of other bonds, whereas the bonds in the latter deforms
independently. The SB-PD was later extended into two types, namely the the
ordinary state based peridynamics (OSB-PD) and the non-ordinary state based
peridynamics (NOSB-PD). In all the above types of PD, horizons sizes are
commonly required to be constant to avoid spurious wave reflections and ghost
forces between particles. However, in many applications, the spatial
distribution of the particles with changing horizon sizes is necessary, e.g.
adaptive refinement, multiscale modelling and multibody analysis. In other
words, in order to achieve acceptable accuracy, the entire numerical model has
to be discretised with respect to the highest particle resolution locally
required, and the smallest horizon size used accordingly. This is
computationally expensive and undesirable. The restriction of horizons being
position independent practically reduces the efficiency of PD.

In this paper, we aim to remove the issue of varying horizons and ghost force
by developing a new PD formulation. The new approach is based on the concept
of horizon and dual-horizon. Though peridynamics has been developed for
different types of physical fields, e.g. thermal field and fluid field, we
confine the present work to solving solid mechanics problems. The content of
the paper is outlined as follows. §2 begins with stating the phenomenon of the
ghost forces. In §3, the horizon and dual-horizon are introduced. New motion
equations with varying horizons are derived based on the dual-horizon concept.
The balance of linear momentum and the balance of angular momentum of the
present PD formulation are proved. The applications of the new formulation for
BB-PD, OSB-PD and NOSB-PD are described in details in §4. In §5, three
numerical examples are presented to validate the present method.

##  2 Ghost force and spurious wave reflection in peridynamics

In the original PD theory, the forces exerted on a particle is the summation
of all the pairwise forces from the particles falling inside the horizon of
that particle. When the horizon sizes are set constant for all particles, the
bond force is always pairwise. On the other hand, it should be noted that the
size of the particles that represent the mass quantity can vary [2, 15].
Spurious wave reflections emerge when the horizon sizes vary. We illustrate
the mechanism of the phenomenon as follows. Consider a particle
$\mathbf{x}^{\prime}$ falling inside the horizon of particle $\mathbf{x}$, see
Fig.1. Let $\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}$ denote the force
vector acting on particle $\mathbf{x}$ due to particle $\mathbf{x}^{\prime}$,
where the first subscript $\mathbf{x}$ indicates $\mathbf{x}$ being the object
of force and the second subscript $\mathbf{x}^{\prime}$ indicates
$\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}$ being the source of force from
$\mathbf{x}^{\prime}$. Take the Fig. 1 for an example, as
$\mathbf{x}^{\prime}\in H_{\mathbf{x}}$, the force
$\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}\neq 0$, $\mathbf{x}\notin
H_{\mathbf{x}^{\prime}}$ due to unequal size of horizons, hence
$\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}=0$. When computing the reactive
forces for particle $\mathbf{x}$, the bond forces
$\mathbf{F}_{\mathbf{x}\mathbf{x}^{\prime}}=\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}-\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}=-\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}$,
which is added to $\mathbf{F}_{\mathbf{x}}$. Likewise, when computing forces
for particle $\mathbf{x}^{\prime}$, particle $\mathbf{x}$ exerts no force on
$\mathbf{x}^{\prime}$ as $\mathbf{x}$ is not inside the horizon of
$\mathbf{x}^{\prime}$. Consequently, the bond force
$\mathbf{F}_{\mathbf{x}\mathbf{x}^{\prime}}$ only exists unilaterally, which
is known as the “ghost force” [1] resulting in an unbalanced internal force.
The balance of linear momentum and balance of angular momentum in this case
are violated, and hence yields spurious wave reflections in PD simulations.

Figure 1: Force vector in peridynamics with varying horizons. Reactive force
$-\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}$ on $\mathbf{x}$ due to
$\mathbf{x}$ exerted on $\mathbf{x}^{\prime}$ as $\mathbf{x}^{\prime}\in
H_{\mathbf{x}}$.

Efforts have been made by researchers to explore the possibility of making PD
suitable for nonuniform spatial discretisation. Bobaru et. al studied the
convergence and adaptive refinement in 1D peridynamics [16] and multi-scale
modeling in 2D BB-PD [17]. Dipasquale et al. recently [18] introduced a
trigger based on the damage state of the material for 2D refinement of BB-PD.
Refined PD was developed in [19, 16, 17, 18, 20], however it is restricted to
BB-PD formulation and to the authors knowledge, has only been used for one-
and two-dimensional problems. In these works, the spurious wave reflection or
ghost force problem is not solved, and the refinement is performed by checking
the spurious reflection is within an acceptable range compared to the
magnitude of the whole wave. Recently, the partial stress was proposed in [1]
to remove the ghost force for varying horizons. However, the method imposes
certain restrictions to the deformation of the body and requires the
computation of partial stress which is complicated. It to certain extent
impairs the simplicity of the original PD, especially for BB-PD and OSB-PD.
Besides, the method does not completely remove the ghost force but with a
small residual which is believed to be acceptable. Though the slice method was
devised in the same work, it is applicable only to piecewise constant sizes of
horizons, and requires additional computation to enforce the consistency
between particles close to the interface.

##  3 Governing equations based on horizon and dual-horizon

In this section, the concept of horizon and dual-horizon will be introduced to
formulate the balance equations for particles with varying horizons. It is
applied to all peridynamic formulations.

###  3.1 Horizon and dual-horizon

We will begin by restating the original concept of horizon in peridynamics. In
peridynamics theory, the particles interact with each other within a finite
distance. A particle is considered to have an influence over other particles
within a small domain centering that particle. The radius of the domain is
known as “horizon”, see Fig.2. Particle $\mathbf{x}$ (thick solid line) is
included in the horizons of particles $\mathbf{x}_{1}$, $\mathbf{x}_{2}$,
$\mathbf{x}_{3}$ and $\mathbf{x}_{4}$ (thin solid line), however is not
included in the horizons of particles $\mathbf{x}_{5}$ and $\mathbf{x}_{6}$
(dashed line). The concept of horizon can be compared to the concept of “nodal
support” in meshless methods.

Figure 2: The schematic diagram for horizon and dual-horizon, all circles
above are horizons. The green points
$\\{\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{3},\mathbf{x}_{4}\\}\in
H_{\mathbf{x}}^{\prime}$,whose horizons denote by thin solid line; the red
points $\\{\mathbf{x}_{5},\mathbf{x}_{6}\\}\notin H_{\mathbf{x}}^{\prime}$
,whose horizons denote by dashed line

Horizon
The horizon $H_{\mathbf{x}}$ is the domain where any particle falling inside
will receive the forces exerted by $\mathbf{x}$. Hence, $\mathbf{x}$ will
undertake all the reactive or passive forces from particles in
$H_{\mathbf{x}}$. Hence, the horizon can be viewed as passive force horizon.
The reactive force acting on $\mathbf{x}$ follows Newton’s third law. As shown
in Fig.1, the reactive force $-\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}$
exerted on $\mathbf{x}$ by other particles is in opposite direction of the
forces applied by $\mathbf{x}$ to other particles.

Dual-horizon
Dual-horizon is defined as a union of points whose horizons include
$\mathbf{x}$, denoted by
$H_{\mathbf{x}}^{\prime}=\\{\mathbf{x}^{\prime}:\mathbf{x}\in
H_{\mathbf{x}^{\prime}}\\}$. In the notation of dual-horizon
$H_{\mathbf{x}}^{\prime}$, the superscript prime indicates “dual”, and the
subscript $\mathbf{x}$ denotes the object particle. It can be understood as
the set of the horizons that belong to the particles who can “see”
$\mathbf{x}$ in their horizons. As shown in Fig.2, the dual-horizon with
respect to $\mathbf{x}$ is the union of $\mathbf{x}_{1}$, $\mathbf{x}_{2}$,
$\mathbf{x}_{3}$ and $\mathbf{x}_{4}$, whose horizons are denoted by thin
solid circles. Particles $\mathbf{x}_{5}$ and $\mathbf{x}_{6}$ are not
included in the dual-horizon of $\mathbf{x}$ since their horizons do not
include $\mathbf{x}$. In this case, $\mathbf{x}$ becomes the object “observed”
by the other particles. If $\mathbf{x}$ is within the horizon of
$\mathbf{x}^{\prime}$, then $\mathbf{x}^{\prime}$ has an active or direct
effect on $\mathbf{x}$, corresponding to the passive effect in horizon defined
previously. Particle $\mathbf{x}$ receives the active forces from other
particles in the dual-horizon of $\mathbf{x}$, and in this sense it is
considered as “dual” corresponding to horizon. For any point $\mathbf{x}$, the
shape of $H_{\mathbf{x}}^{\prime}$ is arbitrary, and depends on the sizes and
shapes of horizons as well as the locations of the particles. Note that the
horizon can take other shapes other than circles or spheres.

The bond inside the dual-horizon ($H_{\mathbf{x}}^{\prime}$) is termed as
“dual-bond”, and this is corresponding to the bond in the horizon
($H_{\mathbf{x}}$). It can be seen that the bond of one particle can become
the dual-bond for another particle interacting with it. Note that for each
particle, the bond and the dual-bond are independent from each other; the same
applies to the horizon and dual-horizon. It means the bond and dual-bond can
break independently in the fracture models. For discretisations with varying
horizons, often one particle is within the horizons of other particles but not
vice versa. In this case, a single horizon is not sufficient to define the
interactions between particles. The concept of two horizons proposed here can
solve this issue with a simple and direct physical meaning. It naturally takes
into account the interactions between particles of varying horizon sizes. For
models with constant horizons, horizon and dual-horizon will degenerate to the
horizon in original peridynamics.

Reaction force by horizon and dual-horizon
In our dual-horizon peridynamic formulation, the horizons are differentiated
between how a particle receives and exerts forces with other particles. Under
this new concept, computing the force
$\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}$ between a pair of particles,
denoted as particles $\mathbf{x}$ and $\mathbf{x}^{\prime}$, is determined by
whether $\mathbf{x}$ is within the horizon of $\mathbf{x}^{\prime}$, and
whether $\mathbf{x}^{\prime}$ is inside the dual-horizon of $\mathbf{x}$, and
vice versa for $\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}$. It can be easily
seen that the force density vector has the following property,

$\displaystyle\text{if}\quad\mathbf{x}\in
H_{\mathbf{x}^{\prime}}\,\text{or}\,\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{{}^{\prime}},\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}\neq 0\
,$ $\displaystyle\text{else}\quad\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}=0\
.$ (2)

For any bond between two particles $\mathbf{x}$ and $\mathbf{x}^{\prime}$
belonging to a domain denoted as $\mathbf{\Omega}$ in Fig.3, the direct force
$\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}$ acting on $\mathbf{x}$ due to
$\mathbf{x}^{\prime}$ can be computed by two approaches as follows.
Approach 1 computes the force in terms of $\mathbf{x}$,

$\displaystyle\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}=\begin{cases}\neq
0&\mbox{if }\mathbf{x}^{\prime}\in H_{\mathbf{x}}^{\prime}\\\ 0&\mbox{if
}\mathbf{x}^{\prime}\notin H_{\mathbf{x}}^{\prime}\end{cases}\ ,$

and Approach 2 is formulated with respect to $\mathbf{x}^{\prime}$,

$\displaystyle\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}=\begin{cases}\neq
0&\mbox{if }\mathbf{x}\in H_{\mathbf{x}^{\prime}}\\\
0&\mbox{if}\mathbf{x}\notin H_{\mathbf{x}^{\prime}}\end{cases}\ .$

Figure 3: Double summation of force

For any domain under consideration e.g. the shaded area in Fig.3, the
computation of the direct forces that take place between particles (no
reactive force considered yet) shall undertake all forces from any particle
that belongs to $\mathbf{\Omega}$. There are two approaches to achieve that.
The first approach is by summing all the forces the particles undertake in the
dual-horizon,

$\displaystyle\sum_{\mathbf{x}\in\mathbf{\Omega}}\left(\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}}\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}\,\triangle
V_{\mathbf{x}^{\prime}}\right)\triangle V_{\mathbf{x}}.$ (3)

And the second is to add up all the forces of each particle that it applies to
other particles,

$\displaystyle\sum_{\mathbf{x}^{\prime}\in\mathbf{\Omega}}\left(\sum_{\mathbf{x}\in
H_{\mathbf{x}^{\prime}}}\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}\,\triangle
V_{\mathbf{x}}\right)\triangle V_{\mathbf{x}^{\prime}}\ .$ (4)

Since the total force for any $\mathbf{\Omega}$ is independent of the approach
chosen to compute it, Eq. (3) and (4) shall be equal

$\displaystyle\sum_{\mathbf{x}\in\mathbf{\Omega}}\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}}\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}\,\triangle
V_{\mathbf{x}^{\prime}}\triangle
V_{\mathbf{x}}=\sum_{\mathbf{x}^{\prime}\in\mathbf{\Omega}}\sum_{\mathbf{x}\in
H_{\mathbf{x}^{\prime}}}\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}\,\triangle
V_{\mathbf{x}}\triangle V_{\mathbf{x}^{\prime}}\ .$ (5)

Eq. (5) indicates that changing the summation order from
$\mathbf{x}\to\mathbf{x}^{\prime}$ shall be done together with changing from
$H_{\mathbf{x}^{\prime}}\to H_{\mathbf{x}}^{\prime}$ so as to keep all the
variables consistent. When the discretisation is sufficiently fine, the
summation shall approximate to the integral and Eq. (5) becomes

$\displaystyle\int_{\mathbf{x}\in\mathbf{\Omega}}\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}}\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}\,\mathbf{d}V_{\mathbf{x}^{\prime}}\mathbf{d}V_{\mathbf{x}}=\int_{\mathbf{x}^{\prime}\in\mathbf{\Omega}}\int_{\mathbf{x}\in
H_{\mathbf{x}^{\prime}}}\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}\,\mathbf{d}V_{\mathbf{x}}\mathbf{d}V_{\mathbf{x}^{\prime}}\
.$ (6)

###  3.2 Motion equation for peridynamics with horizon-variable

In the following derivation $\mathbf{x}^{\prime}-\mathbf{x}$ will be denoted
as $\boldsymbol{\xi}$, and
$\mathbf{u}(\mathbf{x}^{\prime},t)-\mathbf{u}(\mathbf{x},t)$ as
$\boldsymbol{\eta}$, hence $\boldsymbol{\xi}+\boldsymbol{\eta}$ represents the
current relative position vector between the particles. The internal forces
that are exerted at each particle from the other particles should include two
parts, namely the forces from the horizon and the forces from the dual-
horizon. The other forces applied to a particle include the body force and the
inertia force. Let $\triangle V_{\mathbf{x}}$ denote the volume associates to
$\mathbf{x}$. The body force for particle $\mathbf{x}$ can be expressed as
$\mathbf{b}(\mathbf{x},t)\triangle V_{\mathbf{x}}$, where
$\mathbf{b}(\mathbf{x},t)$ is the body force density. The inertia is denoted
by $\rho\mathbf{\ddot{u}}(\mathbf{x},t)\triangle V_{\mathbf{x}}$, where $\rho$
is the density associated to $\mathbf{x}$. At any time $t$, for
$\mathbf{x}^{\prime}$ in dual-horizon of $\mathbf{x}$, the force vector of
$\tilde{\mathbf{f}}_{\mathbf{x}\mathbf{x}^{\prime}}$ is defined as

$\displaystyle\tilde{\mathbf{f}}_{\mathbf{x}\mathbf{x}^{\prime}}:=\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}(\boldsymbol{\eta},\boldsymbol{\xi})\cdot\triangle
V_{\mathbf{x}}\cdot\triangle V_{\mathbf{x}^{\prime}}\ ,$ (7)

where
$\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}(\boldsymbol{\eta},\boldsymbol{\xi})$
is the force density in the traditional peridynamics with unit of force per
volume squared; $\tilde{\mathbf{f}}_{\mathbf{x}\mathbf{x}^{\prime}}$ is the
force vector acting on particle $\mathbf{x}$ due to the attraction or
repulsion from $\mathbf{x}^{\prime}$. The forces from the dual-horizon are
active forces as they are applied to $\mathbf{x}$. The total force applied to
$\mathbf{x}$ from its dual-horizon, $H_{\mathbf{x}}^{\prime}$, can be computed
by,

$\displaystyle\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}}\tilde{\mathbf{f}}_{\mathbf{x}\mathbf{x}^{\prime}}=\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}}\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}(\boldsymbol{\eta},\boldsymbol{\xi})\cdot\triangle
V_{\mathbf{x}}\cdot\triangle V_{\mathbf{x}^{\prime}}\ .$ (8)

For any $\mathbf{x}^{\prime}$ inside the horizon of $\mathbf{x}$, the force
$\tilde{\mathbf{f}}_{\mathbf{x}^{\prime}\mathbf{x}}$ acting on
$\mathbf{x}^{\prime}$ due to $\mathbf{x}$ is defined as

$\displaystyle\tilde{\mathbf{f}}_{\mathbf{x}^{\prime}\mathbf{x}}:=\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}(-\boldsymbol{\eta},-\boldsymbol{\xi})\cdot\triangle
V_{\mathbf{x}^{\prime}}\cdot\triangle V_{\mathbf{x}}$ (9)

where
$\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}(-\boldsymbol{\eta},-\boldsymbol{\xi})$
is the force density per volume squared;
$\tilde{\mathbf{f}}_{\mathbf{x}^{\prime}\mathbf{x}}$ is the force vector
acting on particle $\mathbf{x}^{\prime}$ due to the attraction or repulsion
from $\mathbf{x}$ ( $\mathbf{x}^{\prime}$ inside horizon of $\mathbf{x}$). The
total force $\mathbf{x}$ exerted on other particles from horizon
$H_{\mathbf{x}}$ is the summation of
$\tilde{\mathbf{f}}_{\mathbf{x}^{\prime}\mathbf{x}}$

$\displaystyle\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\tilde{\mathbf{f}}_{\mathbf{x}^{\prime}\mathbf{x}}=\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}(-\boldsymbol{\eta},-\boldsymbol{\xi})\cdot\triangle
V_{\mathbf{x}}\cdot\triangle V_{\mathbf{x}^{\prime}}$ (10)

Based on Newton’s third law, the total force $\mathbf{x}$ undertaken in
$H_{\mathbf{x}}$ takes opposite direction

$\displaystyle-\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\tilde{\mathbf{f}}_{\mathbf{x}^{\prime}\mathbf{x}}=-\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}(-\boldsymbol{\eta},-\boldsymbol{\xi})\cdot\triangle
V_{\mathbf{x}}\cdot\triangle V_{\mathbf{x}^{\prime}}$ (11)

By summing over all forces on particle $\mathbf{x}$, including inertia force,
body force, active force in Eq. (8) and passive force in Eq. (11), we obtain
the equations of motion

$\displaystyle\rho\mathbf{\ddot{u}}(\mathbf{x},t)\triangle
V_{\mathbf{x}}=\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}}\tilde{\mathbf{f}}_{\mathbf{x}\mathbf{x}^{\prime}}+\left(-\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\tilde{\mathbf{f}}_{\mathbf{x}^{\prime}\mathbf{x}}\right)+\mathbf{b}(\mathbf{x},t)\triangle
V_{\mathbf{x}}\ .$ (12)

Substituting Eqs. (8) and (11) into Eq. (12) leads to

$\displaystyle\rho\mathbf{\ddot{u}}(\mathbf{x},t)\triangle
V_{\mathbf{x}}=\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}}\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}(\boldsymbol{\eta},\boldsymbol{\xi})\triangle
V_{\mathbf{x}^{\prime}}\triangle V_{\mathbf{x}}-\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}(-\boldsymbol{\eta},-\boldsymbol{\xi})\triangle
V_{\mathbf{x}^{\prime}}\triangle
V_{\mathbf{x}}+\mathbf{b}(\mathbf{x},t)\triangle V_{\mathbf{x}}\ .$ (13)

As the volume $\triangle V_{\mathbf{x}}$ associated to particle $\mathbf{x}$
is independent of the summation, we can eliminate $\triangle V_{\mathbf{x}}$
in Eq. (13), yielding the governing equation based on $\mathbf{x}$:

$\displaystyle\rho\mathbf{\ddot{u}}(\mathbf{x},t)=\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}}\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}(\boldsymbol{\eta},\boldsymbol{\xi})\triangle
V_{\mathbf{x}^{\prime}}-\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}(-\boldsymbol{\eta},-\boldsymbol{\xi})\triangle
V_{\mathbf{x}^{\prime}}+\mathbf{b}(\mathbf{x},t)\ .$ (14)

When the discretisation is sufficiently fine, the summation is approximating
to the integration of the force on the dual-horizon and horizon,

$\displaystyle\lim_{\triangle V_{\mathbf{x}^{\prime}}\to
0}\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}}\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}(\boldsymbol{\eta},\boldsymbol{\xi})\triangle
V_{\mathbf{x}^{\prime}}=\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}}\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}(\boldsymbol{\eta},\boldsymbol{\xi})\,\mathbf{d}V_{\mathbf{x}^{\prime}}$
(15)

and

$\displaystyle\lim_{\triangle V_{\mathbf{x}^{\prime}}\to
0}\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}(-\boldsymbol{\eta},-\boldsymbol{\xi})\triangle
V_{\mathbf{x}^{\prime}}=\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}(-\boldsymbol{\eta},-\boldsymbol{\xi})\,\mathbf{d}V_{\mathbf{x}^{\prime}}\
.$ (16)

Thus the integration form of the equation of motion in peridynamics with dual
horizon is given as

$\displaystyle\rho\mathbf{\ddot{u}}(\mathbf{x},t)=\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}}\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}(\boldsymbol{\eta},\boldsymbol{\xi})\,\mathbf{d}V_{\mathbf{x}^{\prime}}-\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}(-\boldsymbol{\eta},-\boldsymbol{\xi})\,\mathbf{d}V_{\mathbf{x}^{\prime}}+\mathbf{b}(\mathbf{x},t)\
.$ (17)

Eq. (17) is similar to Eq. (1) of the original peridynamics theory. When the
horizons are set constant, i.e. both horizon and the dual-horizon are equal,
the integrations in Eq. (17) degenerate to the original peridynamics theory.
It means the traditional peridynamics can be viewed as a special case of the
present dual horizon peridynamics.

About the implementation of the present peridynamic formulation, for any
particle $\mathbf{x}$, the force density
$\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}(\boldsymbol{\eta},\boldsymbol{\xi})$
in $H_{\mathbf{x}}^{\prime}$ can be determined when calculating the force in
$H_{\mathbf{x}^{\prime}}$ for particle $\mathbf{x}^{\prime}$. Therefore, it is
not necessary to know exactly the dual-horizon geometry and the formulation
can be implemented with minor modification of any peridynamic codes.

###  3.3 Proof of basic physical principles

####  3.3.1 Balance of linear momentum

The internal forces shall satisfy the balance of linear momentum for any
bounded body $\mathbf{\Omega}$ given by

$\displaystyle\int_{\mathbf{\Omega}}(\rho\mathbf{\ddot{u}}(\mathbf{x},t)-\mathbf{b}(\mathbf{x},t))\mathbf{d}V_{\mathbf{x}}$
$\displaystyle=\int_{\mathbf{\Omega}}\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}}\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}(\boldsymbol{\eta},\boldsymbol{\xi})\,\mathbf{d}V_{\mathbf{x}^{\prime}}\mathbf{d}V_{\mathbf{x}}-\int_{\mathbf{\Omega}}\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}(-\boldsymbol{\eta},-\boldsymbol{\xi})\,\mathbf{d}V_{\mathbf{x}^{\prime}}\mathbf{d}V_{\mathbf{x}}$
$\displaystyle=\mathbf{0}\ .$ (18)

For simplicities, let $\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}$ represent
$\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}(\boldsymbol{\eta},\boldsymbol{\xi})$,
and $\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}$ the
$\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}(-\boldsymbol{\eta},-\boldsymbol{\xi})$.

Proof: For convenience, we start with the discrete form of the Eq. (14) as

$\displaystyle\sum_{\mathbf{x}\in\mathbf{\Omega}}(\rho\mathbf{\ddot{u}}(\mathbf{x},t)-\mathbf{b}(\mathbf{x},t))\triangle
V_{\mathbf{x}}$
$\displaystyle=\sum_{\mathbf{x}\in\mathbf{\Omega}}\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}}\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}\,\triangle
V_{\mathbf{x}^{\prime}}\triangle
V_{\mathbf{x}}-\sum_{\mathbf{x}\in\mathbf{\Omega}}\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}\,\triangle
V_{\mathbf{x}^{\prime}}\triangle V_{\mathbf{x}}$ (19)

The first term on RHS of Eq. (19) can be substituted with Eq. (5) by changing
the summation domain and therefore we will get,

$\displaystyle\sum_{\mathbf{x}\in\mathbf{\Omega}}(\rho\mathbf{\ddot{u}}(\mathbf{x},t)-\mathbf{b}(\mathbf{x},t))\triangle
V_{\mathbf{x}}=\sum_{\mathbf{x}^{\prime}\in\mathbf{\Omega}}\sum_{\mathbf{x}\in
H_{\mathbf{x}^{\prime}}}\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}\,\triangle
V_{\mathbf{x}^{\prime}}\triangle
V_{\mathbf{x}}-\sum_{\mathbf{x}\in\mathbf{\Omega}}\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}\,\triangle
V_{\mathbf{x}^{\prime}}\triangle V_{\mathbf{x}}\ ,$ (20)

The above transformation is based on the idea that
$\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}$ acting on $\mathbf{x}$ is
computed by using dual-horizon ($H_{\mathbf{x}}^{\prime}$) of $\mathbf{x}$.
And $\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}$ is calculated by using
horizon ($H_{\mathbf{x}^{\prime}}$) of $\mathbf{x}^{\prime}$. Note that in
both calculations, $\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}$ remains
unchanged and only the definition domain changes. Since the dummy variables
can be relabeled for $\mathbf{x}\leftrightarrow\mathbf{x}^{\prime}$ in the
first term on RHS, thus the first term has the same expression as the second
term, i.e.

$\displaystyle\sum_{\mathbf{x}\in\mathbf{\Omega}}(\rho\mathbf{\ddot{u}}(\mathbf{x},t)-\mathbf{b}(\mathbf{x},t))\triangle
V_{\mathbf{x}}=\sum_{\mathbf{x}\in\mathbf{\Omega}}\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}\,\triangle
V_{\mathbf{x}}\triangle
V_{\mathbf{x}^{\prime}}-\sum_{\mathbf{x}\in\mathbf{\Omega}}\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}\,\triangle
V_{\mathbf{x}^{\prime}}\triangle V_{\mathbf{x}}=0$

As the forces are differentiated here between active and passive force, i.e.
an active force from $\mathbf{x}$ corresponds to a passive force of
$\mathbf{x}^{\prime}$, the forces is always pairwise but with opposite
direction and of the same magnitude. Hence, it is natural that the balance of
linear momentum is satisfied.

####  3.3.2 Balance of angular momentum

Let $\mathbf{y}$ be the deformation vector state field defined by

$\displaystyle\mathbf{y}[\mathbf{x},t]\langle\boldsymbol{\xi}\rangle=\mathbf{y}(\mathbf{x}^{\prime},t)-\mathbf{y}(\mathbf{x},t)\qquad\forall\mathbf{x}\in\mathbf{\Omega},\boldsymbol{\xi}=\mathbf{x}^{\prime}-\mathbf{x},\mathbf{x}^{\prime}\in
H_{\mathbf{x}},t\geq 0\,$ (21)
$\displaystyle\mathbf{u}(\mathbf{x},t)=\mathbf{y}(\mathbf{x},t)-\mathbf{x}\ ,$
(22)

where $\mathbf{y}(\mathbf{x},t)$ refers to the current configuration
coordinate for $\mathbf{x}$ in material configuration,
$\mathbf{u}(\mathbf{x},t)$ is the displacement for $\mathbf{x}$,
$H_{\mathbf{x}}$ is the horizon. Thus
$\mathbf{y}\langle\mathbf{x}^{\prime}-\mathbf{x}\rangle$ is the image of the
bond $\mathbf{x}^{\prime}-\mathbf{x}$ under the deformation.

To satisfy the balance of angular momentum for any bounded body
$\mathbf{\Omega}$,it is required that:

$\displaystyle\int_{\mathbf{\Omega}}$
$\displaystyle\mathbf{y}\times(\rho\mathbf{\ddot{u}}(\mathbf{x},t)-\mathbf{b}(\mathbf{x},t))\mathbf{d}V_{\mathbf{x}}$
$\displaystyle=\int_{\mathbf{\Omega}}\mathbf{y}\times\left(\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}}\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}(\boldsymbol{\eta},\boldsymbol{\xi})\,\mathbf{d}V_{\mathbf{x}^{\prime}}-\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}(-\boldsymbol{\eta},-\boldsymbol{\xi})\,\mathbf{d}V_{\mathbf{x}^{\prime}}\right)\mathbf{d}V_{\mathbf{x}}$
$\displaystyle=\int_{\mathbf{\Omega}}\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}}\mathbf{y}\times\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}(\boldsymbol{\eta},\boldsymbol{\xi})\,\mathbf{d}V_{\mathbf{x}^{\prime}}\mathbf{d}V_{\mathbf{x}}-\int_{\mathbf{\Omega}}\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\mathbf{y}\times\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}(-\boldsymbol{\eta},-\boldsymbol{\xi})\,\mathbf{d}V_{\mathbf{x}^{\prime}}\mathbf{d}V_{\mathbf{x}}$
$\displaystyle=\mathbf{0}$ (23)

For simplicities, let $\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}$ represent
$\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}(\boldsymbol{\eta},\boldsymbol{\xi})$,
and $\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}$ present
$\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}(-\boldsymbol{\eta},-\boldsymbol{\xi})$.

Proposition: In the dual-horizon peridynamics, suppose a constitutive model of
the form

$\displaystyle\mathbf{f}=\hat{\mathbf{f}}(\mathbf{y},\Lambda)$ (24)

where $\hat{\mathbf{f}}:\mathcal{V}\to\mathcal{V}\ $ is bounded and Riemann-
integrable on $H$ and $\mathcal{V}$ is the vector state; $\Lambda$ denotes all
variables other than the current deformation vector state that may depend on
for some particular material.

If

$\displaystyle\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\mathbf{y}\langle\mathbf{x}^{\prime}-\mathbf{x}\rangle\times\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}\,\mathbf{d}V_{\mathbf{x}^{\prime}}=0\qquad\forall\mathbf{y}\in\mathcal{V}\
,$ (25)

or in discrete form

$\displaystyle\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\mathbf{y}\langle\mathbf{x}^{\prime}-\mathbf{x}\rangle\times\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}\,\triangle
V_{\mathbf{x}^{\prime}}=0\qquad\forall\mathbf{y}\in\mathcal{V}\ ,$ (26)

where
$\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}=\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}(\mathbf{u}(\mathbf{x}^{\prime},t)-\mathbf{u}(\mathbf{x},t),\mathbf{x}^{\prime}-\mathbf{x})$
, the force vector density acting on $\mathbf{x}^{\prime}$, and
$H_{\mathbf{x}}$ is the horizon of $\mathbf{x}$. Then, the balance of angular
momentum, Eq. (23) holds for any deformation of $\mathbf{\Omega}$ for any
given constitutive model.

Proof: As the summation over the domain $\mathbf{\Omega}$ is equivalent to the
integrand of the same expression over that domain, we use the discrete form as

$\displaystyle\sum_{\mathbf{x}\in\mathbf{\Omega}}\mathbf{y}\times(\rho\mathbf{\ddot{u}}(\mathbf{x},t)-\mathbf{b}(\mathbf{x},t))\triangle
V_{\mathbf{x}}$
$\displaystyle=\sum_{\mathbf{x}\in\mathbf{\Omega}}(\mathbf{x}+\mathbf{u})\times(\rho\mathbf{\ddot{u}}(\mathbf{x},t)-\mathbf{b}(\mathbf{x},t))\triangle
V_{\mathbf{x}}$
$\displaystyle=\sum_{\mathbf{x}\in\mathbf{\Omega}}\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}}(\mathbf{x}+\mathbf{u})\times\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}\,\triangle
V_{\mathbf{x}^{\prime}}\triangle
V_{\mathbf{x}}-\sum_{\mathbf{x}\in\mathbf{\Omega}}\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}(\mathbf{x}+\mathbf{u})\times\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}\,\triangle
V_{\mathbf{x}^{\prime}}\triangle V_{\mathbf{x}}$ (27)

The first term on the RHS of Eq. (3.3.2) can be rewritten according to Eq. (5)
by changing the summation domain and thus it becomes

$\displaystyle\sum_{\mathbf{x}\in\mathbf{\Omega}}\mathbf{y}\times(\rho\mathbf{\ddot{u}}(\mathbf{x},t)-\mathbf{b}(\mathbf{x},t))\triangle
V_{\mathbf{x}}$
$\displaystyle=\sum_{\mathbf{x}^{\prime}\in\mathbf{\Omega}}\sum_{\mathbf{x}\in
H_{\mathbf{x}^{\prime}}}(\mathbf{x}+\mathbf{u})\times\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}\,\triangle
V_{\mathbf{x}^{\prime}}\triangle
V_{\mathbf{x}}-\sum_{\mathbf{x}\in\mathbf{\Omega}}\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}(\mathbf{x}+\mathbf{u})\times\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}\,\triangle
V_{\mathbf{x}^{\prime}}\triangle V_{\mathbf{x}}\ .$ (28)

Swapping the dummy variables for the first term on RHS between
$\mathbf{x}\leftrightarrow\mathbf{x}^{\prime}$ in the double summation, then
the first term takes the same summation to the second term on RHS

$\displaystyle\sum_{\mathbf{x}\in\mathbf{\Omega}}\mathbf{y}\times(\rho\mathbf{\ddot{u}}(\mathbf{x},t)-\mathbf{b}(\mathbf{x},t))\triangle
V_{\mathbf{x}}$
$\displaystyle=\sum_{\mathbf{x}\in\mathbf{\Omega}}\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}(\mathbf{x}^{\prime}+\mathbf{u}^{\prime})\times\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}\,\triangle
V_{\mathbf{x}}\triangle
V_{\mathbf{x}^{\prime}}-\sum_{\mathbf{x}\in\mathbf{\Omega}}\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}(\mathbf{x}+\mathbf{u})\times\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}\,\triangle
V_{\mathbf{x}^{\prime}}\triangle V_{\mathbf{x}}$
$\displaystyle=\sum_{\mathbf{x}\in\mathbf{\Omega}}(\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}((\mathbf{x}^{\prime}+\mathbf{u}^{\prime})-(\mathbf{x}+\mathbf{u}))\times\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}\,\triangle
V_{\mathbf{x}^{\prime}})\triangle V_{\mathbf{x}}$
$\displaystyle=\sum_{\mathbf{x}\in\mathbf{\Omega}}(\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}(\mathbf{y}^{\prime}-\mathbf{y})\times\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}\,\triangle
V_{\mathbf{x}^{\prime}})\triangle V_{\mathbf{x}}$
$\displaystyle=\sum_{\mathbf{x}\in\mathbf{\Omega}}(\sum_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\mathbf{y}\langle\mathbf{x}^{\prime}-\mathbf{x}\rangle\times\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}\,\triangle
V_{\mathbf{x}^{\prime}})\triangle V_{\mathbf{x}}=\mathbf{0}\,$ (29)

or in an integration form for sufficiently fine discretisation as

$\displaystyle\int_{\mathbf{x}\in\mathbf{\Omega}}\mathbf{y}\times(\rho\mathbf{\ddot{u}}(\mathbf{x},t)-\mathbf{b}(\mathbf{x},t))\mathbf{d}V_{\mathbf{x}}=\mathbf{0}\
.$ (30)

Therefore, the angular momentum over the entire analysis domain is satisfied.
$\square$

Figure 4: force vector of BB-PD,OSB-PD and NOSB-PD

It can be seen that the conservation of angular somehow depends only on the
horizon; the dual-horizon is not involved. The latter is only needed in the
Eq. (14) and Eq. (17). The bond-based, ordinary based and non-ordinary based
peridynamics all satisfy the angular momentum in the original horizon concept
(see proof in [15]). This conclusion can be also illustrated for the BB-PD and
OS-PD since the internal forces $\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}$
and $\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}$ are parallel to the bond
vector in the current configuration, see Fig.4.

##  4 Dual-horizon peridynamics

The dual-horizon formulation will now be applied to all existing peridynamics,
namely BB-PD, OSB-PD and NOSB-PD. The calibration of constitutive parameters
with respect to the continuum model and some issues concerning the
implementations will be discussed for 2D and 3D problems.

###  4.1 Dual-horizon bond based peridynamics

In the bond-based peridynamics theory [4, 8], the pair force function is
expressed as

$\displaystyle\mathbf{f}(\boldsymbol{\eta},\boldsymbol{\xi})=\frac{\partial
w(\boldsymbol{\eta},\boldsymbol{\xi})}{\partial\boldsymbol{\eta}}\qquad\forall\boldsymbol{\eta},\boldsymbol{\xi}$
(31)

For a micro-elastic homogenous and isotropic material,
$\mathbf{f}(\boldsymbol{\eta},\boldsymbol{\xi})$ can be specified as

$\displaystyle\mathbf{f}(\boldsymbol{\eta},\boldsymbol{\xi})=cs\cdot\frac{\boldsymbol{\eta}+\boldsymbol{\xi}}{\|{\boldsymbol{\eta}+\boldsymbol{\xi}}\|}\
,$ (32)

where $c$ is the micro-modulus, and $s$ is the bond stretch calculated by

$\displaystyle
s=\frac{\|\boldsymbol{\eta}+\boldsymbol{\xi}\|-\|\boldsymbol{\xi}\|}{\|\boldsymbol{\xi}\|}\
.$ (33)

The energy per unit volume in the body at a given time $t$ is given by [8]

$\displaystyle
W=\frac{1}{2}\int_{H_{\mathbf{x}}}w(\boldsymbol{\eta},\boldsymbol{\xi})\mathbf{d}V_{\boldsymbol{\xi}}\
.$ (34)

For the BB-PD theory, by enforcing the strain energy density being equal to
the strain energy density in the classical theory of elasticity [4, 8], and by
setting the influence function $\omega(\|{\boldsymbol{\xi}}\|)=1$, we obtain

$\displaystyle C(\delta)$ $\displaystyle=\frac{3E}{\pi\delta^{3}(1-\nu)}\quad$
$\displaystyle\mathrm{plane\,stress}$ $\displaystyle C(\delta)$
$\displaystyle=\frac{3E}{\pi\delta^{3}(1+\nu)(1-2\nu)}\quad$
$\displaystyle\mathrm{plane\,strain}$ $\displaystyle C(\delta)$
$\displaystyle=\frac{3E}{\pi\delta^{4}(1-2\nu)}\quad$
$\displaystyle\mathrm{3D}\ .$ (35)

Note that the value of $C(\delta)$ takes half of the micro-modulus $c$ used in
the horizon-constant BB-PD since the bond energy for varying horizon is
determined by both horizon and dual-horizon. It can be easily seen that when
the horizon takes the same value as the dual-horizon, the bond energy between
two particles is reduced to that of the original BB-PD.

Let $w_{0}(\boldsymbol{\xi})=C(\delta)s_{0}^{2}(\delta)\xi/2$ denote the work
required to break a single bond, where $s_{0}(\delta)$ is the critical bond
stretch. By breaking half of all the bonds connected to a given particle along
the fracture surface and equalizing the breaking bonds energy with the energy
release rate $G_{0}$ [18, 4], we can get the expression between the energy
release rate $G_{0}$ and the critical bond stretch $s_{0}(\delta)$:

$\displaystyle s_{0}(\delta)$ $\displaystyle=\sqrt{\frac{4\pi
G_{0}}{9E\delta}}\quad$ $\displaystyle\mathrm{plane\,stress}$ $\displaystyle
s_{0}(\delta)$ $\displaystyle=\sqrt{\frac{5\pi G_{0}}{12E\delta}}\quad$
$\displaystyle\mathrm{plane\,strain}$ $\displaystyle s_{0}(\delta)$
$\displaystyle=\sqrt{\frac{5G_{0}}{6E\delta}}\quad$ $\displaystyle\mathrm{3D}\
.$ (36)

Both the micro-modulus $C(\delta)$ and the critical stretch $s_{0}(\delta)$
are derived from the local continuum mechanics theory, and they depend on the
horizon radii for variable horizons.

In the implementation of the bond-based peridynamics, fracture is introduced
by removing particles from the neighbour list once the bond stretch exceeds
the critical bond stretch $s_{0}$. In order to specify whether a bond is
broken or not, a history-dependent scalar valued function $\mu$ is introduced
[8],

$\displaystyle\mu(t,\boldsymbol{\xi})=\begin{cases}1\,\qquad\mathrm{if}\,s(t^{\prime},\boldsymbol{\xi})<s_{0}\,\,\text{for
all }\,0\leq t^{\prime}\leq t,\\\ 0\,\qquad\mathrm{otherwise.}\end{cases}$
(37)

The local damage at $\mathbf{x}$ is defined as

$\displaystyle\phi(\mathbf{x},t)=1-\frac{\int_{H_{\mathbf{x}}}\mu(\mathbf{x},t,\boldsymbol{\xi})\mathbf{d}V_{\xi}}{\int_{H_{\mathbf{x}}}\mathbf{d}V_{\xi}}\
.$ (38)

The damage formulation of Eq. (38) is also applicable to OSB-PD. For any
particle with dual-horizon ($H_{\mathbf{x}}^{\prime}$) and horizon
($H_{\mathbf{x}}$), the active force
$\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}$ and the passive force
$\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}$ in Eq. (14) or (17) are computed
by the following expressions, respectively

$\displaystyle\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}=C(\delta_{\mathbf{x}^{\prime}})\cdot
s_{\mathbf{x}\mathbf{x}^{\prime}}\cdot\frac{\boldsymbol{\eta}+\boldsymbol{\xi}}{\|{\boldsymbol{\eta}+\boldsymbol{\xi}}\|},\qquad\forall\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}$ (39)
$\displaystyle\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}=C(\delta_{\mathbf{x}})\cdot
s_{\mathbf{x}\mathbf{x}^{\prime}}\cdot\frac{-(\boldsymbol{\eta}+\boldsymbol{\xi})}{\|{\boldsymbol{\eta}+\boldsymbol{\xi}}\|},\qquad\forall\mathbf{x}^{\prime}\in
H_{\mathbf{x}}\ ,$ (40)

where $C(\delta_{\mathbf{x}^{\prime}})$ and $C(\delta_{\mathbf{x}})$ are the
micro-modulus based on $\delta_{\mathbf{x}^{\prime}}$ and
$\delta_{\mathbf{x}}$ computed from Eq. (4.1) respectively, and
$s_{\mathbf{x}\mathbf{x}^{\prime}}$ is the stretch between particles
$\mathbf{x}$ and $\mathbf{x}^{\prime}$.

###  4.2 Dual-horizon ordinary state based peridynamics

The concept of “state” for peridynamics was firstly introduced in [15]. A
state of order $m$ is a function

$\underline{\mathbf{A}}:H\to\mathbf{T}_{m},\boldsymbol{\xi}\mapsto\underline{\mathbf{A}}\langle\boldsymbol{\xi}\rangle.$
(41)

where $H$ is the horizon domain, $\mathbf{T}_{m}$ denotes the set of all
tensors of order $m$, $\langle\square\rangle$ indicate the vector to which a
state operates. In this paper, the scalar state and vector state are use if
not otherwise specified.

The 2D OSB-PD formulation of bond force can be derived by following the
similar procedure as that was described in 3D. Both 2D and 3D OSB-PD can be
written in a unified expression taking into account the dimension number as

$\displaystyle\underline{t}\langle\boldsymbol{\xi}\rangle=\frac{nK\theta}{m}\omega\langle\boldsymbol{\xi}\rangle\cdot\xi+\frac{n(n+2)G}{m}\omega\langle\boldsymbol{\xi}\rangle\underline{e}^{d}\langle\boldsymbol{\xi}\rangle\
,$ (42)

where $n\in\\{2,3\\}$ is the dimensional number, $K$ is the bulk modulus, $G$
the shear modulus, $\displaystyle\xi=\|\boldsymbol{\xi}\|$,
$\underline{e}^{d}\langle\boldsymbol{\xi}\rangle=\underline{e}\langle\boldsymbol{\xi}\rangle-\displaystyle{\frac{\theta\xi}{n}}$,
$m=\int_{H}\omega\langle\boldsymbol{\xi}\rangle\,\boldsymbol{\xi}\cdot\boldsymbol{\xi}\,\mathbf{d}V_{\boldsymbol{\xi}}$,
$\displaystyle{\theta=\frac{n}{m}}\int_{H}{\omega\langle\boldsymbol{\xi}\rangle\,\xi\cdot\underline{e}\langle\boldsymbol{\xi}\rangle\mathbf{d}V_{\boldsymbol{\xi}}}$,
$\underline{e}\langle\boldsymbol{\xi}\rangle=\|\boldsymbol{\xi}+\boldsymbol{\eta}\|-\|\boldsymbol{\xi}\|$
and $\mathbf{d}V_{\xi}$ are the area in 2D, and volume in 3D, respectively.
For any particle $\mathbf{x}$, $\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}$
and $\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}$ in Eqs. (14) and (17) can be
calculated by

$\displaystyle\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}=\underline{t}\langle\boldsymbol{\xi}\rangle\cdot\frac{\boldsymbol{\eta}+\boldsymbol{\xi}}{\|{\boldsymbol{\eta}+\boldsymbol{\xi}}\|},\qquad\forall\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}\ ,$ (43)

and

$\displaystyle\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}=\underline{t}^{\prime}\langle-\boldsymbol{\xi}\rangle\cdot\frac{-\boldsymbol{\eta}-\boldsymbol{\xi}}{\|{\boldsymbol{\eta}+\boldsymbol{\xi}}\|},\qquad\forall\mathbf{x}^{\prime}\in
H_{\mathbf{x}}\ .$ (44)

Here $\boldsymbol{\xi}=\mathbf{x}^{\prime}-\mathbf{x}$, and $\mathbf{x}$ and
$\mathbf{x}^{\prime}$ are the coordinate vectors for $\mathbf{x}$ and particle
$\mathbf{x}^{\prime}$ in the material configuration, respectively.

###  4.3 Dual-horizon non-ordinary state based peridynamics

NOSB-PD uses the deformation gradient from classic continuum mechanics. It
offers the possibility to consistently include existing constitutive models.
With the use of the shape tensor proposed in the NOSB-PD, the fundamental
tensors in continuum mechanics can be easily introduced in peridynamics,
including the deformation gradient tensor or velocity gradient [21, 22]. The
shape tensor for particle $\mathbf{x}$ is defined as

$\displaystyle\mathbf{K}_{\mathbf{x}}$
$\displaystyle=\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\omega(\boldsymbol{\xi})\cdot\boldsymbol{\xi}\otimes\boldsymbol{\xi}\mathbf{d}V_{\mathbf{x}^{\prime}}\,$
(45)

where $\boldsymbol{\xi}=\mathbf{x}^{\prime}-\mathbf{x}$ is the bond in the
reference configuration, $H_{\mathbf{x}}$ is the horizon for particle
$\mathbf{x}$ and $\omega(\boldsymbol{\xi})$ is the influence function. The
deformation tensor for particle $\mathbf{x}$ is given by

$\displaystyle\mathbf{F}_{\mathbf{x}}=\frac{\partial{\mathbf{y}}}{\partial{\mathbf{x}}}=\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\omega(\boldsymbol{\xi})\mathbf{y}\langle\boldsymbol{\xi}\rangle\otimes\boldsymbol{\xi}\;\mathbf{d}V_{\mathbf{x}^{\prime}}\,\cdot\mathbf{K}_{\mathbf{x}}^{-1}\
,$ (46)

where $\mathbf{y}:=\mathbf{y}(\mathbf{x},t)$ are the spatial coordinates,
$\mathbf{x}$ are the material coordinates, and
$\boldsymbol{\xi}=\mathbf{x}^{\prime}-\mathbf{x}$. The spatial velocity
gradient $\mathbf{L}_{\mathbf{x}}$ for particle $\mathbf{x}$ is defined in the
current configuration using the chain rule,

$\displaystyle\mathbf{L}_{\mathbf{x}}:=\frac{\partial{\mathbf{v}}}{\partial{\mathbf{y}}}=\frac{\partial{\mathbf{v}}}{\partial{\mathbf{x}}}\cdot\frac{\partial{\mathbf{x}}}{\partial{\mathbf{y}}}=\dot{\mathbf{F}}_{\mathbf{x}}\mathbf{F}_{\mathbf{x}}^{-1}\
,$ (47)

where
$\mathbf{v}:=\mathbf{v}(\mathbf{x},t)=\displaystyle{\frac{\partial\mathbf{y}}{\partial
t}}(\mathbf{x},t)$ is the velocity vector and $\dot{\mathbf{F}}_{\mathbf{x}}$
is the rate of deformation gradient. The non-local form of the velocity
gradient can be written as

$\displaystyle\mathbf{L}_{\mathbf{x}}:$
$\displaystyle=\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\omega(\boldsymbol{\xi})\,{\mathbf{v}}\langle\boldsymbol{\xi}\rangle\otimes\boldsymbol{\xi}\;\mathbf{d}V_{\mathbf{x}^{\prime}}\,\cdot\mathbf{K}_{\mathbf{x}}^{-1}\cdot\left(\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\omega(\boldsymbol{\xi})\,\mathbf{y}\langle\boldsymbol{\xi}\rangle\otimes\boldsymbol{\xi}\;\mathbf{d}V_{\mathbf{x}^{\prime}}\,\cdot\mathbf{K}_{\mathbf{x}}^{-1}\right)^{-1}$
$\displaystyle=\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\omega(\boldsymbol{\xi})\,\mathbf{v}\langle\boldsymbol{\xi}\rangle\otimes\boldsymbol{\xi}\;\mathbf{d}V_{\mathbf{x}^{\prime}}\,\cdot\mathbf{K}_{\mathbf{x}}^{-1}\cdot\mathbf{K}_{\mathbf{x}}\left(\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\omega(\boldsymbol{\xi})\,\mathbf{y}\langle\boldsymbol{\xi}\rangle\otimes\boldsymbol{\xi}\;\mathbf{d}V_{\mathbf{x}^{\prime}}\right)^{-1}\,\cdot$
$\displaystyle=\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\omega(\boldsymbol{\xi})\,{\mathbf{v}}\langle\boldsymbol{\xi}\rangle\otimes\boldsymbol{\xi}\;\mathbf{d}V_{\mathbf{x}^{\prime}}\cdot\left(\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\omega(\boldsymbol{\xi})\,\mathbf{y}\langle\boldsymbol{\xi}\rangle\otimes\boldsymbol{\xi}\;\mathbf{d}V_{\mathbf{x}^{\prime}}\right)^{-1}$
(48)

where $\boldsymbol{\xi}=\mathbf{x}^{\prime}-\mathbf{x}$,
$\mathbf{v}\langle\boldsymbol{\xi}\rangle=\mathbf{v}(\mathbf{x}^{\prime},t)-\mathbf{v}(\mathbf{x},t)$
and
$\mathbf{y}\langle\boldsymbol{\xi}\rangle=\mathbf{y}(\mathbf{x}^{\prime},t)-\mathbf{y}(\mathbf{x},t)$.
It can been seen that the velocity gradient does not require shape tensor,
which means the shape tensor is not necessary to define velocity gradient. The
incremental deformation gradient can also be obtained without shape tensor
$\mathbf{K}$.

Let $\triangle\mathbf{u}$ denotes the displacement increment with respect to
the previous time step. The incremental spatial deformation gradient
$\mathbf{C}_{\mathbf{x}}$ in continuum mechanics is defined as

$\displaystyle\mathbf{C}_{\mathbf{x}}:=\frac{\partial{(\triangle\mathbf{u})}}{\partial{\mathbf{y}}}=\frac{\partial{(\triangle\mathbf{u})}}{\partial{\mathbf{x}}}\cdot\frac{\partial{\mathbf{x}}}{\partial{\mathbf{y}}}=\frac{\partial{(\triangle\mathbf{u})}}{\partial{\mathbf{x}}}\cdot\mathbf{F}_{\mathbf{x}}^{-1}\
,$ (49)

The non-local counterpart is given by

$\displaystyle\mathbf{C}_{\mathbf{x}}:=\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\omega(\boldsymbol{\xi})\triangle\mathbf{u}\langle\boldsymbol{\xi}\rangle\otimes\boldsymbol{\xi}\;\mathbf{d}V_{\mathbf{x}^{\prime}}\cdot(\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\omega(\boldsymbol{\xi})\mathbf{y}\langle\boldsymbol{\xi}\rangle\otimes\boldsymbol{\xi}\;\mathbf{d}V_{\mathbf{x}^{\prime}})^{-1}\quad\
.$

Many material models in non-ordinary state based peridynamic accounting for
the geometrical and material non-linear problems are solved on the incremental
spatial deformation gradient. The non-ordinary bond force is computed based on
the Piola-Kirchhoff as

$\displaystyle\underline{\mathbf{T}}_{\mathbf{x}^{\prime}\mathbf{x}}=\omega(\boldsymbol{\xi})\,\mathbf{P}_{\mathbf{x}}\cdot\mathbf{K}_{\mathbf{x}}^{-1}\cdot\boldsymbol{\xi}$
(50)
$\displaystyle\mathbf{P}_{\mathbf{x}}=J_{\mathbf{x}}\sigma_{\mathbf{x}}\mathbf{F}_{\mathbf{x}}^{-T},J_{\mathbf{x}}=\det{\mathbf{F}_{\mathbf{x}}}$
(51)

where $\boldsymbol{\xi}=\mathbf{x}^{\prime}-\mathbf{x}$,$\sigma_{\mathbf{x}}$
is Cauchy stress tensor.

Let us define the deformation gradient $\mathbf{F}_{\mathbf{x}}^{\prime}$

$\displaystyle\mathbf{F}_{\mathbf{x}}^{\prime}:=\mathbf{F}_{\mathbf{x}}\cdot\mathbf{K}_{\mathbf{x}}=\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\omega(\boldsymbol{\xi})\mathbf{y}\langle\boldsymbol{\xi}\rangle\otimes\boldsymbol{\xi}\;\mathbf{d}V_{\mathbf{x}^{\prime}}\,\cdot\mathbf{K}_{\mathbf{x}}^{-1}\cdot\mathbf{K}_{\mathbf{x}}=\int_{\mathbf{x}^{\prime}\in
H_{\mathbf{x}}}\omega(\boldsymbol{\xi})\mathbf{y}\langle\boldsymbol{\xi}\rangle\otimes\boldsymbol{\xi}\;\mathbf{d}V_{\mathbf{x}}\
.$ (52)

Therefore

$\displaystyle
J_{\mathbf{x}}=\det{\mathbf{F}_{\mathbf{x}}}=\frac{\det{\mathbf{F}_{\mathbf{x}}^{\prime}}}{\det\mathbf{K}_{\mathbf{x}}}\
.$ (53)

By substituting Eqs. (53), (52) and (51) to (50), we obtain

$\displaystyle\underline{\mathbf{T}}_{\mathbf{x}^{\prime}\mathbf{x}}$
$\displaystyle=\omega(\boldsymbol{\xi})\,J_{\mathbf{x}}\sigma_{\mathbf{x}}\mathbf{F}_{\mathbf{x}}^{-T}\cdot\mathbf{K}_{\mathbf{x}}^{-1}\cdot\boldsymbol{\xi}$
$\displaystyle=\omega(\boldsymbol{\xi})\,\frac{\det{\mathbf{F}_{\mathbf{x}}^{\prime}}}{\det{\mathbf{K}_{\mathbf{x}}}}\cdot\sigma_{\mathbf{x}}\cdot(\mathbf{F}_{\mathbf{x}}^{\prime}\cdot\mathbf{K}_{\mathbf{x}}^{-1})^{-T}\cdot\mathbf{K}_{\mathbf{x}}^{-1}\cdot\boldsymbol{\xi}$
$\displaystyle=\omega(\boldsymbol{\xi})\,\frac{\det{\mathbf{F}_{\mathbf{x}}^{\prime}}}{\det{\mathbf{K}_{\mathbf{x}}}}\cdot\sigma_{\mathbf{x}}\cdot\mathbf{F}_{\mathbf{x}}^{\prime-T}\cdot\mathbf{K}_{\mathbf{x}}\cdot\mathbf{K}_{\mathbf{x}}^{-1}\cdot\boldsymbol{\xi}$
$\displaystyle=\omega(\boldsymbol{\xi})\,\frac{\det{\mathbf{F}_{\mathbf{x}}^{\prime}}}{\det{\mathbf{K}_{\mathbf{x}}}}\cdot\sigma_{\mathbf{x}}\cdot\mathbf{F}_{\mathbf{x}}^{\prime-T}\cdot\boldsymbol{\xi}\
.$ (54)

The shape tensor is only needed for
$J_{\mathbf{x}}=\det{\mathbf{F}_{\mathbf{x}}}$, while the computation of the
bond force density does not require the shape tensor. For any particle with
horizon $H_{\mathbf{x}}$ and dual-horizon $H^{\prime}_{\mathbf{x}}$ , the
active force $\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}$ and the passive
force $\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}$ in motion Eq. (14) or (17)
can be calculated by

$\displaystyle\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}=\underline{\mathbf{T}}_{\mathbf{x}\mathbf{x}^{\prime}}\qquad\forall\mathbf{x}^{\prime}\in
H_{\mathbf{x}}^{\prime}\ ,$ (55)

and

$\displaystyle\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}=\underline{\mathbf{T}}_{\mathbf{x}^{\prime}\mathbf{x}}\qquad\forall\mathbf{x}^{\prime}\in
H_{\mathbf{x}}\ ,$ (56)

where $\boldsymbol{\xi}=\mathbf{x}^{\prime}-\mathbf{x}$, and $\mathbf{x}$ and
$\mathbf{x}^{\prime}$ are the coordinate vectors for particle $\mathbf{x}$ and
particle $\mathbf{x}^{\prime}$ in the material configuration, respectively.

The summary of three kinds of dual-horizon peridynamics see Table.1.

| original Peridynamics | dual-horizon Peridynamics
---|---|---
Model | $\int_{H_{\mathbf{x}}}(\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}-\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}})\mathbf{d}V_{\mathbf{x}^{\prime}}$ | $\int_{H_{\mathbf{x}}^{\prime}}\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}\mathbf{d}V_{\mathbf{x}^{\prime}}-\int_{H_{\mathbf{x}}}\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}\mathbf{d}V_{\mathbf{x}^{\prime}}$
BB-PD |  $\displaystyle\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}-\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}=c\cdot s_{\mathbf{x}\mathbf{x}^{\prime}}\cdot\mathbf{n},\forall\mathbf{x}^{\prime}\in H_{\mathbf{x}}$ $\displaystyle\mbox{where }\mathbf{n}=\frac{\boldsymbol{\eta}+\boldsymbol{\xi}}{\|{\boldsymbol{\eta}+\boldsymbol{\xi}}\|}$ |  $\displaystyle\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}=C(\delta_{\mathbf{x}^{\prime}})\cdot s_{\mathbf{x}\mathbf{x}^{\prime}}\cdot\mathbf{n},\forall\mathbf{x}^{\prime}\in H_{\mathbf{x}}^{\prime};$ $\displaystyle\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}=C(\delta_{\mathbf{x}})\cdot s_{\mathbf{x}\mathbf{x}^{\prime}}\cdot(-\mathbf{n}),\forall\mathbf{x}^{\prime}\in H_{\mathbf{x}}$ $\displaystyle\mbox{where }\mathbf{n}=\frac{\boldsymbol{\eta}+\boldsymbol{\xi}}{\|{\boldsymbol{\eta}+\boldsymbol{\xi}}\|}$
OSB-PD |  $\displaystyle\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}=\underline{t}\langle\boldsymbol{\xi}\rangle\cdot\mathbf{n},$ $\displaystyle\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}=\underline{t}^{\prime}\langle-\boldsymbol{\xi}\rangle\cdot(-\mathbf{n}),\forall\mathbf{x}^{\prime}\in H_{\mathbf{x}}$ $\displaystyle\mbox{where }\mathbf{n}=\frac{\boldsymbol{\eta}+\boldsymbol{\xi}}{\|{\boldsymbol{\eta}+\boldsymbol{\xi}}\|}$ |  $\displaystyle\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}=\underline{t}\langle\boldsymbol{\xi}\rangle\cdot\mathbf{n},\forall\mathbf{x}^{\prime}\in H_{\mathbf{x}}^{\prime};$ $\displaystyle\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}=\underline{t}^{\prime}\langle-\boldsymbol{\xi}\rangle\cdot(-\mathbf{n}),\forall\mathbf{x}^{\prime}\in H_{\mathbf{x}}$ $\displaystyle\mbox{where }\mathbf{n}=\frac{\boldsymbol{\eta}+\boldsymbol{\xi}}{\|{\boldsymbol{\eta}+\boldsymbol{\xi}}\|}$
NOSB-PD |  $\displaystyle\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}=\underline{\mathbf{T}}_{\mathbf{x}\mathbf{x}^{\prime}},$ $\displaystyle\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}=\underline{\mathbf{T}}_{\mathbf{x}^{\prime}\mathbf{x}},\forall\mathbf{x}^{\prime}\in H_{\mathbf{x}}$ |  $\displaystyle\mathbf{f}_{\mathbf{x}\mathbf{x}^{\prime}}=\underline{\mathbf{T}}_{\mathbf{x}\mathbf{x}^{\prime}},\forall\mathbf{x}^{\prime}\in H_{\mathbf{x}}^{\prime};$ $\displaystyle\mathbf{f}_{\mathbf{x}^{\prime}\mathbf{x}}=\underline{\mathbf{T}}_{\mathbf{x}^{\prime}\mathbf{x}},\forall\mathbf{x}^{\prime}\in H_{\mathbf{x}}$
Table 1: Comparison of the original Peridynamics and dual-horizon Peridynamics

##  5 Numerical Examples

###  5.1 Two-dimensional wave reflection in a rectangular plate

Consider a rectangular plate with dimensions of 0.1$\times$0.04 $\mbox{m}^{2}$
(see Fig.5). The Young’s modulus, density and the Poisson’s ratio used for the
plate are $E=1$, $\rho=1$ and $\nu=0$, respectively.Note that this is 1-D
model. The initial displacement applied to the plate is described by the
following equation

$\displaystyle u_{0}(x,y)=0.0002\exp[-(\frac{x}{0.01})^{2}],\quad
v_{0}(x,y)=0,\quad x\in[0,0.1],y\in[-0.02,0.02]\ ,$ (57)

where $u_{0}$ and $v_{0}$ denote the displacement in the $x$ and $y$
directions respectively. The wave speed is $v=\sqrt{E/\rho}=1\,\mbox{m/s}$. At
any time step, the $L_{2}$ error in the displacement is given by

$\displaystyle{\|\text{err}\|}_{L_{2}}=\frac{\|\mathbf{u}^{h}-\mathbf{u}_{\text{analytic}}\|}{\|\mathbf{u}_{\text{analytic}}\|}\
,$ (58)

with

$\displaystyle\|\mathbf{u}\|=\left(\int_{\Omega_{0}}\mathbf{u}\cdot\mathbf{u}\,d\Omega_{0}\right)^{\frac{1}{2}}\
.$

Figure 5: setup of the plate

Three models for solving this problem, namely model A, B and C, are
devised,see Figs.6(a) and 6(b). In model A, the plate is discretised with 4000
particles and all particles have the same horizon sizes of 0.001 m, see
Fig.6(a). For models B and C, we adopted a discretisation of varying horizons
as shown in Fig.6(b). Verlet integration method is adopted.The horizon radius
associated to each particle is set as 3 times the particle sizes for all
models. Therefore, along the interface between the coarse and fine
discretization, the horizon sizes vary. The minimal particle size of 3 models
is $\triangle x=$5e-4m, yielding a critical time increment $\triangle
t_{max}=\triangle x/v$=5e-4s. The physical computational time step is set 0.5
second with the time increment $\triangle t_{max}$=5e-4s.

As BB-PD can not simulate the model with Poisson ratio $\nu=0$, the OSB-PD
model is adopted. Model A was solved with the OSB-PD with constant horizons.
Model B is the original OSB-PD with variable horizon (without additional
treatment for ghost force), and model C is our dual-horizon formulation for
OSB-PD with variable horizon. The parameters used for the three models are
listed in Table.2.

Model | $\triangle x=\triangle y$ | $10^{3}\cdot\delta$ | Particle numbers
---|---|---|---
A | 0.001 | 3 | 4000
B | 0.001,0.0005 | 3/1.5 | 8628
C | 0.001,0.0005 | 3/1.5 | 8628
Table 2: Values of the Peridynamic parameters for the mixed model

(a) Particles distribution of model A

(b) Particles distribution of model B and C

Figure 6: The discretizations of model A,B and C.

Fig.7,8,9 show the displacement in the x-direction along the x-coordinate of
all particles. The red star points are particles with $\Delta
x=0.0005\mbox{m}$, the blue dot particles refer to $\Delta x=0.001\mbox{m}$.

Figure 7: Displacement wave $u_{x}$ versus $x$ coordinate of model A at step
650 Figure 8: Displacement wave $u_{x}$ versus $x$ coordinate of model B at
step 650 Figure 9: Displacement wave $u_{x}$ versus $x$ coordinate of model C
at step 650

The displacement in model C (present dual-horizon PD) is almost identical to
that of model A (PD with constant horizon),as shown in Fig.7 and Fig.9.
Spurious wave reflections are observed for Model B, as shown in Fig.8. At step
650, several wave peaks are observed in Model B and the maximum displacement
is lower than that of Models A and C. Therefore, the displacement wave in
Model B were affected by spurious wave reflection and the results deviate
greatly from results of A and C.

We compute the $L_{2}$ errors of the wave profile along the line-$CD$ (see
Fig. 5). The $L_{2}$ error (see Table. 3) of the present peridynamic
formulation, can achieve the almost the same accuracy as the conventional
peridynamics with constant horizon.

Model | step 100 | step 200 | step 650 | step 990
---|---|---|---|---
A | 0.0410 | 0.0452 | 0.1277 | 0.2498
B | 0.1149 | 0.1738 | 0.7638 | 1.3233
C | 0.0447 | 0.0500 | 0.1253 | 0.2871
Table 3: $L_{2}$ error of displacement for the initial condition (Gauss
distribution of displacement) before and after the wave reflection

Displacement and velocity of monitor points
Two monitor points (A,B) are selected to evaluate the spurious wave reflection
of Models B and C as shown in Fig. 5. The displacement and velocity curves
show the spurious wave reflection exists in model B while not exists in model
C, see Fig.10.

(a) $x$ displacement of model C

(b) $y$ displacement of model C

(c) $x$ velocity of model C

(d) $y$ velocity of model C

(e) $x$ displacement of model B

(f) $y$ displacement of model B

(g) $x$ velocity of model B

(h) $y$ velocity of model B

Figure 10: The displacement and velocity curve of monitor points $A$ and $B$
for model B and C

###  5.2 Kalthof-Winkler experiment

To study the performance of the present formulation for fracture modeling, the
Kalthof-Winkler experiment is exploited. The geometry and model used for the
test is depicted in Fig.11, the thickness of the specimen is set as 0.01m. In
the experiment, the evolution of crack pattern was observed to be dependent on
the impact loading velocity. For a plate made of steel 18Ni1900 subjected to
an impact loading at the speed of $v_{0}$=32 m/s, brittle fracture was
observed [23]. The crack propagates from the end of the initial crack at an
angle around 70∘ vs. the initial crack direction. In [23], the BB-PD with
constant horizon was used to test this example. For convenience of comparison,
the present dual-horizon formulation is also applied to the bond-based
peridynamics to test the example. The material parameters used are the same as
in [23], i.e. the elastic modulus
$E=190\,\mbox{GPa}$,$\rho=7800\,\mbox{kg/m}^{3}$,$\nu=0.25$ and the energy
release rate $G_{0}=6.9e4\,\mbox{J/m}^{2}$. The impact loading was imposed by
applying an initial velocity at $v_{0}=22\,\mbox{m/s}$ to the first three
layers of particles in the domain as shown in Fig.11;all other boundaries are
free. The plate is discretized with two different particle sizes, namely
$\triangle x_{\text{coarse}}=1.5625\mbox{e-}3\,\mbox{m}$ for the coarse
subdomain and $\triangle x_{\text{dense}}=0.5\triangle
x_{\text{coarse}}=7.8125\mbox{e-}4\,\mbox{m}$ for the fine subdomain located
in the left down corner of the model, see Fig.11. Along the thickness of the
plate (in direction perpendicular to the plane surface), four layers of
particles in the coarse subdomain and eight layers in the fine subdomain are
employed. The total number of particles is 57968.

Figure 11: Kalthof-Winkler’s experimental setup

The crack propagation speed is computed by

$\displaystyle
V_{l-0.5}=\frac{\|\mathbf{x}_{l}-\mathbf{x}_{l-1}\|}{t_{l}-t_{l-1}}\ ,$ (59)

where $\mathbf{x}_{l}$ and $\mathbf{x}_{l-1}$ are the positions of the crack
tip at the times $t_{l}$ and $t_{l-1}$ respectively. An expression for the
Rayleigh wave speed $c_{R}$ [24] is given as

$\displaystyle\frac{c_{R}}{c_{s}}=\frac{0.87+1.12\nu}{1+\nu}\ ,$ (60)

where $\nu$ is the Poisson’s ration, $c_{s}=\sqrt{\mu/\rho}$ is the shear wave
speed and $\mu$ is the shear modulus. The crack starts to propagate at 26.3
$\mu$s. The highest crack speed reached is 1530 m/s, about 54.4% of the
Rayleigh speed(2799.2 m/s). The average crack speed is 1077 m/s for the fine
subdomain and 1094 m/s for the coarse subdomain. During the crack propagation,
the crack paths in the fine and coarse subdomains are nearly symmetrical, as
shown in Figs. 12, 13 and 14. It can be seen from Fig.15 that the crack speed
is very close to that predicted by the PD formulation using constant horizons
($\triangle x_{\text{uniform}}=1.5625e\mbox{-3\,m}$). The crack propagation
initiated at an angle of 65.7∘ in the fine subdomain with respect to the
original crack and 65.8∘ in fine subdomain. Therefore, it can be concluded
that with the present formulation, the transition from fine to coarse of
horizons has almost no influence on the crack propagation speed.

Figure 12: The crack pattern of Kalthof-Winkler plate by the present dual
horizon PD at step 350 Figure 13: The crack pattern of Kalthof-Winkler
simulation by the present dual horizon PD at step 650 Figure 14: The crack
pattern of Kalthof-Winkler simulation by the present dual horizon PD at step
875 Figure 15: The crack speed of Kalthof-Winkler simulation by the present
dual horizon PD

###  5.3 Adaptively refined peridynamics

The present dual-horizon Peridynamic formulation provides the possibility for
adaptive-refinement within a simple and unified framework. Two examples of
adaptively refined peridynamics will be tested in this section: the Kalthof
Winkler experiment in section 5.2 and plate with pre-crack subjected to
traction. The threshold for the adaptive refinement is determined by both the
damage-state criteria and energy state, as proposed in [18]. The adaptive
refinement procedure consists of two steps.
(1) Search for the particles that exceed the threshhold values. Note that the
refinement is not only applied to the particles above the threshold but also
to the neighbouring particles. In this way, it can be ensure that the crack
tip always remains inside the refined zone.
(2) Split the particle, named as parent particle, into small particles, named
as child particles. The properties that will be mapped from the parent
particles to the child particles include the mass, volume, coordinate,
displacement and velocity.
The method to split particles for structured discretization is shown in Fig.
16. Splitting the particle results in halving the maximum time step. In all
examples, we restrict each particle being allowed to split only once in the
entire analysis.

Figure 16: Particle splitting for Adaptive refined peridynamics

####  5.3.1 3D adaptively refined Kalthof-Winkler simulation

We now repeat the Kalthof Winkler experiment from section 5.2. However, the
initial discretization is based on uniform horizon size. The initial particle
size is $\triangle x_{\text{initial}}=$1.5625e-3 m. After refinement, the
particle size is reduced to $\triangle x_{\text{refined}}=$7.8125e-4 m around
the crack. The total particles is increased from 32,768 to 54,860 at the end
of the simulation. A uniform refinement would result 32,768$\times$8=262,144
partilces.

The crack patterns at certain steps are plotted in Figs.17, 18 and 19. The
crack tip is always contained inside the refined zone. According to Fig. 20,
the maximum crack propagation speed reaches 1268 m/s and the average speed is
1047.6 m/s. The crack propagation takes place at an angle of approximately
66.5∘ with respect to the initial crack direction.

Figure 17: The crack pattern of adaptive refined Kalthof-Winkler experiment at
step 420 Figure 18: The crack pattern of adaptive refined Kalthof-Winkler
experiment at step 600 Figure 19: The crack pattern of adaptive refined
Kalthof-Winkler experiment at step 925 Figure 20: The crack speed of adaptive
refined Kalthof-Winkler experiment

####  5.3.2 Plate with pre-crack subjected to traction

In the last example, we will test the capability of the new formulation in
modelling the crack branching. Hence, a plate with pre-crack subject to
traction is considered as shown in Fig.21. The traction remains constant
during the entire simulation. The same example has been well studied in [25]
with BB-PD and also in [26] using the XFEM.

The material parameters of the plate (soda-lime glass) are $E=72\,\mbox{GPa}$,
$\rho=2440\,\mbox{kg/m}^{3}$ and energy release rate being
$G_{0}=135\,\mbox{J/m}^{2}$ [25]. Two models, namely Case 1 and Case 2, are
set up for comparison. Case 1 is solved by using the traditional BB-PD with
constant horizons, while Case 2 is modeled by an adaptively refined BB-PD
using the present dual-horizon formulation. The plate is assumed under plane
stress condition and all simulations are carried out in 2D. The particle sizes
chosen are $5\times 10^{-4}$ m for Case 1 and $1\times 10^{-3}$ m for Case 2,
respectively. Case 2 uses an adaptively refined model, and once a particle is
refined, the minimal particle size will become the same as in Case 1. The
particle number are 16,000 for Case 1 and 4000-6424 for Case 2, respectively.
The Rayleigh speed for the material is 3102 m/s according to Eq. (60). The max
crack speed is 1881.4 m/s (60.6% of $c_{R}$) for Case 1 and 2184.1 m/s (70.4%
of $c_{R}$) for Case 2, respectively. One possible reason for the different
max crack speeds is the deviation of the mapping method.

For Case 1, the crack starts to propagate at 11.9 $\mu$s, and the first crack
branching point $B1$ occurs at 24.5 $\mu$s at the speed of 1043.6 m/s and the
second crack branch point $B2$ takes place at 40.9 $\mu$s at the speed of
1147.1 m/s. It is observed once branching initiates, the crack propagation
speed decreases. This can be explained that the energy release to create more
fracture surfaces decelerates the propagation speed. Afterward, the crack
speed increases. For Case 2, the initial crack started to propagate at 10.3
$\mu$s, first branching point $B1$ at 19.4 $\mu$s with 1247.6 m/s, and second
branching point $B2$ at 34.2 $\mu$s with 1330.6 m/s. The crack pattern is
similar compared with Case 1.

Figure 21: Setup of the pre-cracked plate under traction load Figure 22: The
crack pattern of pre-cracked plate with uniform horizons Figure 23: The crack
pattern of adaptive refined pre-cracked plate Figure 24: The crack speed of
pre-cracked plate under traction

##  6 Conclusions

This paper contributes to the development of peridynamics with varying
horizons. Therefore, the interactions between particles are based on two
independent horizons, passive force from the horizon and the active force from
the dual-horizon. The spurious wave reflection and ghost force in the
conventional peridynamics is naturally eliminated. Based on the new concepts
of horizons, the motion equation of dual-horizon peridynamics is derived. We
show that the balance of linear momentum and angular momentum is satisfied.
The key difference of the present peridynamics formulation from the
conventional one is the way of computing reactions forces. Hence, the dual-
horizon peridynamics can be implemented in any existing peridynamics code with
minimal efforts. Three numerical examples show the present method is free from
spurious wave reflection, and the accuracy is retained along the interface
where horizon sizes undergo sudden changes. The method also shows its
capability for fracture problems including crack branching. A simple
h-adaptive scheme is proposed to improve crack paths resolution, both for 2D
and 3D case. We also intend to develop efficient error estimate for adaptivity
to drive the adaptive refinement. The present method also shows the potential
in multiscale analysis where the models at different length scales can be
bridged by using different horizon settings.

## Acknowledgements

The authors acknowledge the supports from FP7 Marie Curie Actions ITN-INSIST
and IIF-HYDROFRAC (623667), the National Basic Research Program of China (973
Program: 2011CB013800) and NSFC (51474157), the Ministry of Science and
Technology of China (Grant No.SLDRCE14-B-28, SLDRCE14-B-31). The authors would
like to express the gratitude to Dr. Jafar Dashlejeh and Dr. Erkan Oterkus for
the beneficial discussion about Peridynamic theory.

## References

## References

  * [1] S.A. Silling, D.J. Littlewood, and P. Seleson.  Variable horizon in a peridynamic medium.  Technical report, Sandia National Laboratories (SNL-NM), Albuquerque, NM (United States), 2014.
  * [2] S.A. Silling.  Reformulation of elasticity theory for discontinuities and long-range forces.  Journal of the Mechanics and Physics of Solids, 48(1):175–209, 2000\.
  * [3] Xin Lai, Bo Ren, Houfu Fan, Shaofan Li, CT Wu, Richard A Regueiro, and Lisheng Liu.  Peridynamics simulations of geomaterial fragmentation by impulse loads.  International Journal for Numerical and Analytical Methods in Geomechanics, 2015.
  * [4] W. Gerstle, N. Sau, and S. Silling.  Peridynamic modeling of concrete structures.  Nuclear engineering and design, 237(12):1250–1258, 2007.
  * [5] E. Oterkus, E. Madenci, O. Weckner, S. Silling, P. Bogert, and A. Tessler.  Combined finite element and peridynamic analyses for predicting failure in a stiffened composite curved panel with a central slot.  Composite Structures, 94:839–850, 2012.
  * [6] E.T. Moyer and M.J. Miraglia.  Peridynamic solutions for Timoshenko beams.  Engineering, 6:Article ID:46262, 2014.
  * [7] J. O Grady and J. Foster.  Peridynamic plates and flat shells: A non-ordinary, state-based model.  International Journal of Solids and Structures, 51:4572–4579, 2014\.
  * [8] S.A. Silling and E. Askari.  A meshfree method based on the peridynamic model of solid mechanics.  Computers & structures, 83(17):1526–1535, 2005.
  * [9] T. Belytschko, N. Moës, S. Usui, and C. Parimi.  Arbitrary discontinuities in finite elements.  International Journal for Numerical Methods in Engineering, 50(4):993–1013, 2001.
  * [10] T. Rabczuk and T. Belytschko.  Cracking particles: A simplified meshfree method for arbitrary evolving cracks.  International Journal for Numerical Methods in Engineering, 61(13):2316–2343, 2004.
  * [11] I. Babuška and J.M. Melenk.  The partition of the unity finite element method.  Technical Report BN-1185, Institute for Physical Science and Technology, University of Maryland, Maryland, 1995.
  * [12] K.W. Cheng and T.P Fries.  Higher-order xfem for curved strong and weak discontinuities.  International Journal for Numerical Methods in Engineering, DOI: 10.1002/nme.2768.
  * [13] A. Gravouil, N. Moes, and T. Belytschko.  Non-planar 3D crack growth by the extended finite element and level sets - part ii: Level set update.  International Journal for Numerical Methods in Engineering, 53:2569–2586, 2002.
  * [14] S.A. Silling and R.B. Lehoucq.  Convergence of peridynamics to classical elasticity theory.  Journal of Elasticity, 93(1):13–37, 2008.
  * [15] Stewart A Silling, M Epton, O Weckner, J Xu, and E Askari.  Peridynamic states and constitutive modeling.  Journal of Elasticity, 88(2):151–184, 2007.
  * [16] F. Bobaru, M. Yang, L.F. Alves, S.A. Silling, E. Askari, and J. Xu.  Convergence, adaptive refinement, and scaling in 1d peridynamics.  International Journal for Numerical Methods in Engineering, 77(6):852–877, 2009.
  * [17] F. Bobaru and Y.D. Ha.  Adaptive refinement and multiscale modeling in 2d peridynamics.  Journal for Multiscale Computational Engineering, 9(635-659), 2011\.
  * [18] D. Dipasquale, M. Zaccariotto, and U. Galvanetto.  Crack propagation with adaptive grid refinement in 2d peridynamics.  International Journal of Fracture, 190(1-2):1–22, 2014.
  * [19] E. Askari, F. Bobaru, R.B. Lehoucq, M.L. Parks, S.A. Silling, and O. Weckner.  Peridynamics for multiscale materials modeling.  In Journal of Physics: Conference Series, volume 125, page 012078\. IOP Publishing, 2008.
  * [20] K. Yu, X.J. Xin, and K.B. Lease.  A new adaptive integration method for the peridynamic theory.  Modelling and Simulation in Materials Science and Engineering, 19(4):045003, 2011.
  * [21] S.A. Silling and R.B. Lehoucq.  Peridynamic theory of solid mechanics.  Advances in Applied Mechanics, 44:73–168, 2010.
  * [22] T.L. Warren, S.A. Silling, A. Askari, O. Weckner, M.A. Epton, and J. Xu.  A non-ordinary state-based peridynamic method to model solid material deformation and fracture.  International Journal of Solids and Structures, 46:1186–1195, 2008\.
  * [23] S.A. Silling.  Peridynamic modeling of the Kalthoff–Winkler experiment.  Submission for the, 2001.
  * [24] K.F. Graff.  Wave motion in elastic solids.  Clarendon Oxford, 1975.
  * [25] Y.D. Ha and F. Bobaru.  Studies of dynamic crack propagation and crack branching with peridynamics.  International Journal of Fracture, 162:229–244, 2010.
  * [26] T. Belytschko, H. Chen, J. Xu, and G. Zi.  Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment.  International Journal for Numerical Methods in Engineering, 58(12):1873–1905, 2003.