The Multi-Horizon Peridynamics

: We present a novel refinement approach in peridynamics (PD). The proposed approach takes advantage of the PD flexibility in choosing the shape of the horizon by introducing multiple domains (with no intersections) to the nodes of the refinement zone. We will show that no ghost forces are needed when changing the horizon sizes in both subdomains. The approach is applied to both bond-based and state-based peridynamics and verified for a simple wave propagation refinement problem illustrating the efficiency of the method.

approach including advantages and disadvantages. Examples are presented in Section 3 before Section 4 concludes the manuscript.

Multi-horizon peridynamics
The PD equations of motion theoretically introduce no limitation on the integration over their horizon [Silling (2000)]. In this article, we are proposing a change of the horizon shapes for nodes in the interaction regions between domains with different horizon sizes. The PD domain can be refined utilizing the "standard" formulation, which requires a fixed spherical horizon with a constant radius. The refinement cost will increase by the neighborhood size, which is exponentially increasing by the refinement factor (the difference between refined and coarse subdomains). Employing different horizon radii will lead to erroneous results or even instabilities due to ghost forces as illustrated in Fig. 1, where a subdomain at a position X has a horizon radius of R and its neighbor at X ' has a horizon radius of r, where R>r and r<|X-X ' |<R; X will include X ' in its neighborhood but will not be included in the X ' neighborhood. Thus, X will have a force toward X ' (i.e., ghost force) without X ' having any force towards X, causing an imbalance in the bond between X and X ' . For the sake of simplicity, let us assume refinement occurs along a straight line, as illustrated in Fig. 2. The proposed approach suggests enforcing the symmetry of the interaction nodes having different horizon sizes by forcing all the nodes to include themselves to all of their neighbor nodes. For instance, x' in Fig. 2(a) will include all volumes of the subdomains in which x' is included (e.g., x in Fig. 2(a)). Thus, the nodes in the refined domain may have more than one horizon. Note that the interaction between the domains is a one to one relation, thus multi-horizons of the refined nodes cannot coincide. Let us call the part of the multi-horizon that has the same radius as those of refined nodes as the natural horizon (H n ), and the rest of the multi-horizon as the interaction horizon (H i ).
where , , and are the density, the body force, and the displacement of the center of x's horizon respectively. Note, that equation 1 is valid for all PD-types, i.e., bond-based (BB-PD), ordinary state-based (OSB-PD), and non-ordinary state-based (NSB-PD) peridynamics. Let us consider the force vector state T. Its relation to the deformed direction vector state M is given as where can be written as a scalar function over the domain for OSB-PD, and in the special case where ≡ 1, Eq.
(2) represents BB-PD and any other form of function for results in NSB-PD. Eq.
(1) within the refinement zone can be written as where s are the force vector state for the interaction horizon and ′ has to be defined in a way that Eq. (1), and Eq.
(3) have one to one equivalent integral functions on their righthand side. In other words, the multi-horizon formulation has to provide the same acceleration as if the whole domain was refined. This requires where − is the difference between the "standard" horizon and the natural horizon, which is a subset of the interaction horizon since it still has no intersections with the natural horizon. Rewriting Eq. (2) for and utilizing the deformed direction vector state , we obtain = , Substituting Eq. (5) into Eq. (6) Since is constant for both BB-PD and OSB-PD, we can rewrite Eq. (6) as following where the only unknown parameter is , which can be easily computed having the interaction horizon and its equivalent horizon in case of refining the whole domain ( − ). Note that α remains time-independent and therefore, it needs to be computed only at the initial configuration. It is also worth mentioning that NSP-PD can also be implemented as a multihorizon method if Eq. (6) is satisfied. If the refined node is located within a finite number of refinement zones, Eq. (7) can be utilized for each of the interaction horizons individually.

Absence of ghost forces
Ghost forces occur due to violation of Newton's law, which is the case for unsymmetrical interactions of particles, which commonly occur for nodes with different horizon sizes. This is also possible for un-symmetric or dissimilar shapes of the horizons. For models with only spherical horizon shapes for all sub-domains, only differences in the horizon radius between two sub-domains can cause ghost forces in the refinement zone. Multihorizons guarantee the existence of the nodes inside the domain with larger horizon radii, which eliminates ghost forces.

Numerical implementation
Implementing the concept of multi-horizons into an existing PD-code requires three changes: 1. At the initial configuration, the nodes in the refinement zone have to be determined. 2. The parameter α(s) for each node located inside refinement zone should be determined. Note that neighbor nodes may share their representative volume with more than one horizon. 3. After computing the bond forces, the forces of the interaction horizons have to be multiplied with α. Implementing the first two changes usually requires the modification of the data structure in which the horizons and their neighbors are stored.

Numerical example
The computation of α can be cumbersome for complex refinement zones, especially in 2D and 3D, as the intersections of the associated volumes of the neighboring nodes may demand complex geometry computations. For the sake of simplicity, we, present a simple 1D example to demonstrate the performance of the multi-horizon approach. Consider a 1D bar of length 100 mm and two fixed ends applied to a velocity wave excitation in the middle of the bar, as illustrated in Fig. 3. The two ends of the bar have a node distance of 0.1 mm. The distance between nodes and the horizon radius in the middle of the bar is four-time bigger than the ends. The velocity wave has an exponential equation of v =   To study the artificial wave reflection added by refinement zone, the response of the bar is recorded whenever the displacement of the mid node arrives at its peak. Fig. 4 presents the first eight returned velocity waves. The artificial wave reflection due to the refinement zone is relatively small. As illustrated in Fig. 5, the maximum error of 15% is observed after the eights wave reflection.  Fig. 4 compare to the non-refined bar (Colors are the same as Fig. 4) Let us consider now a 0.2 mm node distance in the middle and 0.05 mm node distance at the two ends of the bar. The velocity and its error after the eight wave reflection can be found in Figs. 6 and 7, respectively.   Fig. 6 compare to the non-refined bar (Colors are the same as Fig. 6) Finally, we test a pulse excitation on the same bar (see Fig. 3) where the node distances are 0.05, and 0.2 mm at the end and middle, respectively. Fig. 8 illustrates the first eight wave reflections. Although the error is about 40%, the simulation still remains stable. Note that this error is expected as PD is not well suited for capturing such sharp wave shapes.

Conclusions
A simple approach for PD with variable horizons is proposed, which is important for computational efficiency. The approach avoids ghost forces and related artificial wave reflections. Several problems have been simulated utilizing the proposed method. In summary: • The proposed model can be used in the interaction of any finite number of domains with different node distances and horizon sizes. • The method eliminates the existence of any ghost forces.
• The implementation of the method has almost no effect on the computational cost, and only demands changes in the initiation of the domain(s) configuration. On the downside, the implementation of the method in 2D and 3D requires complex geometrical solver. However, such solvers are open-source and can be combined with the proposed scheme.