Wave reflection and cut‐off frequencies in coupled FE‐peridynamic grids

Reflections are typically observed when pulses propagate across interfaces. Accordingly, spurious reflections might occur at the interfaces between different models used to simulate the same medium. Examples of such coupled models include classical continuum descriptions with molecular dynamics or peridynamic (PD) grids. In this work, three different coupling approaches are implemented to couple bond‐based PDs with finite element (FE) solvers for solid mechanics. It is observed that incorporation of an overlapping zone, over which the coupling between FE and PD occurs, can lead to minimization of the reflected energy compared to a standard force coupling at the FE domain/PD grid interface. However, coupling with other existing methodologies, like the addition of ghost particles, achieves comparable accuracy at lower computational cost. Furthermore, the prudent selection of the discretization parameters is of pivotal importance as they control the high frequency cut‐off limit. Mismatch between the cut‐off frequencies of the different descriptions can lead to unrealistic results.


INTRODUCTION
The computational cost associated with the implementation of the peridynamic (PD) theory is an important issue for large scale analyses. 1 One way to circumvent this drawback is through the incorporation of adaptive refinement techniques near areas where higher resolution is required (eg, near loads, boundaries, crack tips) whilst reducing the number of particles away from areas with stress concentrations. Such an approach leads to a more computationally efficient model and has been analyzed by several authors, eg, Bobaru and Ha 2 and Dipasquale et al. 3 A different approach is through coupling classical continuum elasticity (CCM) with the PD theory. It is preferred that CCM is used in areas where failure is not expected to happen, whereas PD are used to model the behaviour at the vicinity of a failure site, like crack tips, voids, inclusions, notches, etc. The solution of the underlying partial differential equation of classical elasticity is typically approximated using solvers based on the finite element (FE) method. 4 The PD equation of motion on the other hand is usually discretized using a grid of material points/particles and the solution is approximated at their location. 5,6 Because the coupling between these two solutions is commonly realised after the discretization of CCM and PD, it is usually described in the literature as coupling FE with PD grids (see, eg, other works [7][8][9][10]. This terminology is also adopted here. Using such a coupling approach, it is envisaged to combine the computational efficiency of FE with the ability of PD to accommodate discontinuous fields and microstructural characteristics within the solution domain. 11 The advantages of coupling a computational efficient local theory with the capabilities that are offered in finer-scale descriptions (eg, molecular dynamics) have already been identified within multiscale applications. 12 One of the challenges in the coupling methodologies arises when dynamic problems are considered. Due to model mismatch, unrealistic reflections occur when waves propagate across the domains simulated using different models. Substantial effort has been devoted to formulating a methodology that is able to effectively communicate the information between the various descriptions, with several approaches having been proposed to achieve coupling between continuum with atomistic/molecular descriptions. A recent review on this subject can be found in the work of Miller and Tadmor. 13 The coupling methods can be divided into two broad classes: (i) coupling achieved at an interface between the two descriptions of the solid and (ii) coupling being realised gradually over a finite zone, often called overlapping zone. 14,15 It is noted that one of the advantages of overlapping zones is that they can reduce the spurious reflections as well as allow the two models to have different levels of discretization refinement. 14,16,17 Different levels of discretization can be also achieved using interface coupling through the "mortar method" where Lagrange multipliers are used to implicitly impose continuity. 18,19 Many methodologies to couple FE with PDs have also been reported. Without being exhaustive, some contributions are reported here. Liu and Hong 8 have formulated a methodology of coupling FE with PDs using coupling forces through interface elements. In the work of Seleson et al, 20 the PD and the continuum models were blended together using a force-based coupling scheme. The fundamental difference in this approach lies in the fact that the coupled scheme was derived from a single framework instead of using two models and "gluing" them together. 20 The main advantage is that ghost forces are eliminated on the interface. In the works of Kilic and Madenci 9 and Madenci and Oterkus, 21 an overlapping zone is used for the coupling of the two models taking advantage of the FE shape functions. In the work of Han and Lubineau, 22 the Arlequin method, introduced by Dhia and Rateau, 23 is used to couple FE with the nonlocal theory presented by Di Paola et al. 24 Additionally, Han and Lubineau 22 mentioned that PDs is a suitable alternative. In subsequent publications, 1,25 the same authors use the morphing strategy, a modification of the Arlequin method, to couple PDs with FE. It is noted that, in this procedure, a damage index is used to adaptively introduce PDs in areas that are susceptible to fracture. A more complete overview of methodologies to couple FE with PDs can also be found in the work of Zaccariotto et al. 7 In most cases, solutions using coupled FE-PD are restricted to static applications. Recently, Zaccariotto et al 7 used the coupling methodology, introduced in the work of Galvanetto et al, 26 to study dynamic problems as well. Furthermore, applications were extended to 3D problems, illustrating the ability of the coupled model to tackle more complex crack front geometries. The advantage in their methodology lies in its simplicity that does not require determination of blending functions for the two models and can be easily extended to commercial software. Exploiting the computational prowess of graphic devices able to accelerate computations, 27 3D simulations of dynamic crack propagations can become achievable. The existence of spurious reflections due to the coupling procedure was evaluated by transmitting pulses in 1D and 2D scenarios and measuring the amplitude of the reflected wave for different PD horizon values.
In this study, a more systematic analysis regarding the appearance of spurious reflections and the underlying parameters that control the accurate transmission of incident pulses across the coupling interface/zone is presented. Three different coupling methodologies are assessed to evaluate the severity of spurious reflections generated. Emphasis is given on the importance of the cut-off frequency associated with the discretized models. In Section 2, the FE and PD models, used to describe the different domains, are defined. Section 3 provides the formulation of the different coupling approaches considered as well as a description of the problem solved. In total, three approaches are presented, ie, (i) a simple FE-PD coupling using a coincident particle and a node at the interface position, (ii) coupling using additional ghost particles, and (iii) an energy coupling of the two models over an overlapping zone. In Section 4, numerical results of the time domain analyses are presented for a 1D bar. The spurious reflections appearing in each approach are compared for different model parameters. Pulses with different shapes and frequency content are also considered. The coupling approaches presented in Section 3 are also used in Section 5 to couple PD with a continuum model, allowing to perform numerical investigations in the frequency domain and compute reflection and transmission coefficients in an infinite 1D bar. The importance of carefully selecting the discretization parameters to control the model's cut-off frequency is highlighted. Based on the findings from Section 5, one of the coupling approaches is extended to 2D problems in Section 6 to simulate pulses propagating in a plate of infinite dimensions. The unrealistic near trapping of pulses in PD regions due to ineffective coupling with FE is demonstrated, indicating that energy released from, eg, crack propagation might fail to be transferred away from the source. Finally, concluding remarks are included in Section 7.

DEFINITIONS
In this section, the approaches adopted to model a solid will be briefly discussed.

Model 1 -elasticity
First, the local theory of elasticity is briefly presented. In all cases that numerical solutions are pursued, the FE method will be adopted. The problem domain is defined as FE . u and F denote the portions of the boundary , where displacements u d are prescribed and surface forces F d are applied, respectively. Throughout this work, bold notation is used to denote vector fields. The problem is thus stated as follows.
Given the initial conditions u(x,0) and .
u (x, 0) and for x ∈ FE , t ∈ [0,T], find u such thaẗ where u is the displacement vector, is the stress tensor, n is the outward unit vector, and f d are the applied body forces. To enable the FE discretization of FE , the weak (or variational) formulation of Equation (1) is used. Then, the aforementioned problem can be rewritten as follows. Given the initial conditions u(x, 0) and udS (2) for all admissible .
u. In addition, it is (u) = 1 2 (∇u + ∇u T ) and E is the fourth-order elastic tensor. Without loss of generality, in the following, we assume that no body forces and tractions are applied. The displacement field u is approximated by classical FE shape functions as follows: where N np is the total number of mesh nodes in FE , U is the displacement vector at the FE nodal positions, and N i (x) are continuous piecewise polynomial base functions. The discrete form of Equation (2) can be written as follows: where F is the force vector and M and K are the mass and stiffness matrices, respectively, defined as follows: where N en is the total number of elements and i is the element domain. In the following sections, 1D and 2D pulse propagation problems are solved. The approximation of the solution in FE is carried out using two-node linear bar elements for the 1D problems and four-node linear rectangular elements for the 2D problems.

Model 2 -bond-based peridynamics
Bond-based PD is a nonlocal theory introduced by Siling. 28 The theory is based on a partial integrodifferential equation of motion 29 that avoids the use of spatial derivatives. As such, it is suitable for applications where discontinuous displacement fields appear (eg, fracture). In the literature, the PD theory has been used to simulate a wide variety of problems. [30][31][32][33] However, its implementation leads to higher computational requirements compared to FE.
According to the PD theory, a solid body is composed of material points of infinitesimal volume dV called PD particles. Each particle interacts with other particles in its vicinity through a predefined force function f. 21,28 The force function contains the material information along with a length parameter, . This length parameters specifies the extent of the vicinity H x , within which a particle interacts with others. The basic aspects of the PD formulation are briefly summarised here for completeness. We denote with PD the domain where PD theory is used. Following the work of Bobaru et al, 34 the PD equation of motion can be written as follows: where is the density, d is the displacement vector, x is the position vector, b is the body force field, and f is the pairwise force function, derivable from a scalar micropotential w, ie, From the definition of bond-based PD, the pairwise forces are assumed to be collinear, equal in magnitude and opposite to satisfy linear and angular momentum, 28 ie, where = x ′ − x is the relative position vector and = d ′ − d is the relative displacement vector.
Using the definition of a microelastic material, 5 the pairwise force function can be written in the form where s is the bond stretch and k( ) is the bond constant. The value of the bond constant can be evaluated by equating the PD strain energy with its counterpart from classical elasticity. 32 The value of k( ) depends only on the relative position of two particles and describes a linear force-stretch (f − s) particle interaction relationship analogous to a spring idealisation. Nevertheless, the PD equation of motion describes a nonlinear system of equations. In the work of Silling, 28 a linearized version of the PD theory was also derived. The linearized force function reduces to where, C( ) is the second-order micromodulus tensor. 32 Throughout this work, the linearized definition of f is used. However, the coupling approaches presented in the following sections are extendible cases where the original formulation of f is used. The integral of Equation (6) can be approximated numerically as a finite summation using the collocation method described in the work of Bobaru et al 32 as follows: where N is the number of subdomains that are included within the horizon of the ith point, V j is the volume of the jth collocation point, and f i, j is the bond force between points i and j. A more elaborate discretization method has been proposed by Kilic and Madenci,9 where the aforementioned discretization can be extracted as a special case. In Equation (11), either the original or the linearized definition of f can be used. The discretized expression of the PD equation of motion requires the implementation of two correction factors. (i) A surface correction to account for cases where the horizon of a particle is limited due to geometrical boundaries; in the literature, this is often referred to as the "skin effect." 35 The simple method presented in the work of Bobaru et al 32 is implemented here but more refined approaches have been suggested in the literature (see, eg, the work of Madenci and Oterkus 21 ). (ii) A volume correction to account for cases where particles are only partially within the horizon of another particle. Here, the volume correction presented in the work of Liu and Honh 36 is employed.
The variation of the PD kinetic energy can be approximated through a finite summation as follows: where N is the number of particles in PD , . d is the variation of the velocity, and m i = V i is the mass associated with the ith point.
The scalar micropotential w in Equation (7) represents the strain energy that is stored in a single bond and the strain energy density of a material point is calculated by summing the energies of the bonds within its horizon as 37 where the factor 1/2 appears because the energy stored in the bond is distributed equally to each particle that connects. Following a similar approach as in the work of Han and Lubineau 22 and approximating the integral with a finite summation, the variation of the total potential energy in the PD domain can be approximated as follows: If no external forces or body forces are applied in the PD domain, the problem can be stated as follows.
Given the initial conditions d(x, 0) and for every admissible . d. Similar derivations can be found in related works. 16,22,32 It is noted that both the original and the linearized version of PD can be extracted from Equation (15) by making the appropriate substitution of f from Equation (9) or Equation (10).

Time integration
Time integration is very important for all transient problems. Originally, PD were formulated to study dynamic crack propagation in materials under impact or blast loading. 38 Since such phenomena are associated with small time scales and require many time discretization steps, explicit methods have been used widely for the time integration of the PD equation of motion. Among the commonly used time marching techniques are the central difference scheme, implemented in other works, 5,36,39 and the velocity-Verlet algorithm used in the work of Liu and Hong 8 and Ha and Bobaru. 40 Furthermore, time marching performed through forward and backward Euler differences was used in the work of Madenci and Oterkus. 21 In this work, time integration is performed implementing the central difference method 41 (ie, Newmark-with = 0 and = 0.5) for both FE and PD. The displacements of the next time step (n + 1) are computed through Subsequently, next step accelerations are found using the next step displacements U n+1 and d n+1 i from the FE and PD equations of motion. Finally, the velocities of the next step are

COUPLING APPROACHES
To identify the fundamental mechanisms and phenomena associated with spurious wave reflections, a simple case of waves propagating in a 1D bar is used ( Figure 1). The bar domain , consists of the subdomains FE and PD , where = FE ∪ PD and coupl = FE ∩ PD , is the domain where FE-PD coupling takes place. In FE , the classical local theory of elasticity is used, and FE is implemented for the numerical approximation. In PD , the PD model is adopted. Both ends of the bar are contained in FE , whereas PD contains the middle part. In this setting, a propagating wave will cross both an FE-PD and a PD-FE interface. Using this approach has the additional benefit of circumventing the difficulty in PD to define boundary conditions. 7,28,42 Test cases of propagating pulses with different shapes are considered to identify the parameters (eg, PD horizon , effect of discretization size in each model) that are crucial. Using different cases of propagating waves, the ability of each coupling methodology to accurately transfer information across the interface of the two models is studied.
A crucial point of this study is the effect of different coupling approaches on the amount of energy reflected. In total, three different coupling approaches are considered.
i. The first approach (coupling case 1 - Figure 2) involves a simple coupling between the two models. The coupling domain coupl consists of a single node that defines the coupling interface. Force equilibrium and continuity requirements are enforced on the interface to achieve information passing between the two models. Following this approach, it is required that an FE node and a PD particle coincide at the interface location. In this study, this is the only case where such a requirement is enforced. ii. The second approach (coupling case 2 - Figure 2) introduces additional PD particles within the FE mesh that will be referred to as "ghost particles" in the subsequent sections. Recently, an approach very similar to this coupling case was presented in the works of Zaccariotto et al 7 and Galvanetto et al. 26 Similar to coupling case 1, coupl consists of a single node but there is no requirement that this node coincides with a PD particle. Of course, such a requirement could be applied. However, this approach will not be pursued further as realisation of this condition on general 2D and 3D interfaces can be difficult. iii. Lastly, an energy coupling method similar to the Arlequin method (coupling case 3 - Figure 2) developed in the work of Dhia and Rateau 23 is implemented to enforce the coupling. In this case, the two models coexist in the overlapping region and coupl is defined over a region.
For all examples referring to a 1D bar, the value of the bond constant k( ) is assumed to be uniform within the PD horizon = n x PD . Its value is defined as 34

Coupling case 1
Although the first coupling case presented here is simplistic, it can be used as a reference benchmark for comparison with the other two approaches. For convenience, examples for all approaches are described using = 2 x but the formulation remains the same for other values of . Assuming that the interface is at node (eg, = 5 in Figure 2), to ensure continuity between the models at the interface, it is required that The interface node is subjected to forces from the FE model on the left side and the PD model on the right side and the motion of node is thus described from where m = m FE + m PD , m FE , and m PD are the mass contributions from the FE and PD model, respectively, and F FE and F PD are the total forces applied from the FE and the PD model, respectively. Combining Equations (4), (11), (19), and (20), the final system of equations for coupling case 1 is derived. If the linearized formulation of PD is used in PD , it is possible to formulate the stiffness matrix of the PD bonds and write the bar response in the form MÜ + KU = 0. For a single FE, the stiffness matrix is given as follows: The stiffness of each PD bond is given as follows: where the volume correction v i, j has been included. The stiffness matrix of the bar near the interface will be computed as follows: Stiffness contributions from both the PD and the FE domain appear at the interface position (node 5). Furthermore, as also stated in the works of Zaccariotto et al, 7,10 when such couplings are enforced, the bandwidth of the final stiffness matrix changes due to the nonlocal nature of PD.

Coupling case 2
In the second coupling case, ghost particles are inserted within the FE domain ( Figure 2). This approach is very similar to the coupling approach presented in the aforementioned works. 7,10 The number of ghost particles is such that the PD horizon of the first particle in the PD model is not interrupted. Continuity between the two models is enforced by interpolating the nodal displacements of the FE model at the positions of the ghost particles using the FE shape functions.
The introduction of ghost particles does not imply that coincident FE nodes and PD particles on the interface. The forces acting on the interface from the PD model are computed by considering the forces of each individual bond that crosses the coupling interface. Denoting with f i, j the force per unit volume squared particle i exerts on particle j, then the motion of the interface node is described by Combining these time equations (4), (11), (24), and (25), the final system of equations for the second coupling case is formed. Again, using the linearized PD formulation, it is possible to construct the stiffness matrix of the system. Using the numbering in Figure 2, the displacements of the ghost particles can be written in terms of the interpolated nodal displacements as follows: Using the definition of the bond stiffness in Equation (22) and combining Equations (25) and (26), the forces acting on the interface node are 27) and the stiffness matrix becomes Comparing the stiffness matrices of the first two coupling approaches, it is obvious that the nonlocality of PD has also been transferred to the interface node, whereas the horizon of the particles near the interface is not interrupted. It is noted that the stiffness matrix in Equation (28) is valid only for when the PD horizon is = 2 x PD and the ghost particles g1 and g2 lie between nodes ( − 2, − 1) and ( − 1, ), respectively. Of course, similar expressions can be derived for other values of the horizon and ghost particle location.

Coupling case 3
In the third case, coupling is not realised at a single point but occurs along an overlapping region coupl , where the two models coexist ( Figure 2). Coupling between the two domains is enforced in a weak sense through the implementation of Lagrange multipliers. The methodology considered is the one used by Aubertin et al, 16 for continuum to discrete coupling. Similar approaches include the Arlequin method 23 and the bridging domain method. 43 In the work of Aubertin et al, 16 FE were coupled with molecular dynamics. Here, the methodology presented in the aforementioned work 16 is adopted with the difference that the PD theory has been used instead of molecular dynamics.
Coupling in the overlapping region is achieved by introducing a coupling operator that will compare the velocities of the two models. Other coupling operators have been proposed in the literature as well. They enforce the coupling by considering the displacement fields (termed L 2 coupling) or the displacement field and its spatial derivatives (termed H 1 coupling). The interested reader can refer to the works of Dhia and Rateau 23 and Guidault and Belytschko 44 and the references therein for more information. Following the work of Aubertin et al, 16 a "mediator space" M is introduced where the velocities are projected, such that the velocity field is now defined only on discrete points that belong to a subset of PD . The velocities are projected on M using a suitable projection operator Π, a scalar product p, and Lagrange multipliers , such that Similar to the case of molecular dynamics, 16 the discrete nature of PD implies that Lagrange multipliers are defined at the location of PD particles in the overlapping region. The final coupled system can be written as follows. Given the initial conditions u(x,0), d (x, 0), find u and d such that for all admissible  between the two descriptions. In many cases, these functions are called blending functions. 14 They can be defined (using the partition of unity property) through Selection of the appropriate function h(x) to distribute the energies is not straightforward and can be case dependent. 7 In the works of Guidault and Belytschko 44 and Bauman et al, 45 it was demonstrated that, when the L 2 coupling is implemented in the Arlequin method, which is very similar to the coupling implemented here, selection of Heaviside functions to distribute the energy can lead to ill-posed problems. To avoid this, continuous functions were selected to define the energy distribution ( Figure 3). Two cases are used for comparison, ie, (i) a linear function with h (x) = After discretization of Equation (30) and using matrix notation, we arrive at where C FE and C PD are the coupling matrices of the FE and PD, respectively, and a FE and a PD are diagonal weighting matrices. The FE velocities . U are projected on the mediator space M and then interpolated on the location of the PD particles in the overlapping region through the FE shape functions. Matrix C FE contains the shape function values that define this interpolation. C PD is a scalar Boolean matrix associating each Lagrange multiplier with a PD particle.
For the time integration of the third coupling approach, the predictor-corrector procedure presented in the work of Aubertin et al 16 is adopted. It has been shown that this procedure can further reduce the developing of spurious reflections. According to this procedure, each model is first solved individually, disregarding the coupling terms, and then the velocities and accelerations are updated in the overlapping domain using the Lagrange multipliers. The procedure proposed is briefly presented here.
First, the next step displacements U n+1 and d n+1 are calculated using Equations (16a) and (16b). Then, a prediction of the next step accelerationsÜ n+1 , * andd n+1, * is made as follows: Using the predicted accelerations, a prediction of the next step . U n+1, * and . d n+1, * velocities is made using Equations (17a) and (17b). Following the work of Aubertin et al, 16 the Lagrange multipliers are computed using the predicted velocities as follows: Using the Lagrange multipliers calculated from Equation (34), the velocities are adjusted as follows: and the accelerations as follows:

Time-domain analyses
Using the three coupling approaches presented, wave propagation test cases are simulated numerically for the 1D bar illustrated in Figure 1. It is assumed that Young's modulus of the bar is E = 200 GPa, the cross-sectional area is A = 10 −6 m 2 , the density is = 7850 kg/m 3 , and the total length is l = 1.2 m. The bar is divided into three equal parts with the middle one being modelled using the PD theory and the other two with FEs. When the third coupling case is used, the PD domain is extended within the FE domain to create the overlapping region. Test cases with different pulse shapes are carried out. The displacement of the node at the free end of the bar is prescribed to simulate a rightward travelling incident wave. Incident pulses with two different shapes are used: a gaussian with The following parameters are selected for the pulses: a 2 = 4 × 10 −6 , a 3 = 1 × 10 −6 , f 1 = 5 × 10 5 Hz, and f 2 = 5.25 × 10 5 Hz. In both cases, the amplitude of the pulse is a 1 = 10 −4 m. The frequency content of each pulse shape is illustrated in Figure 4.
Coupling case 1 is used first to connect FE with nonlinear PD. Both models have identical discretization lengths (dashed line) are illustrated. Even with a simple approach like this, it is possible to achieve information passing between FE-PD and PD-FE interfaces given that a very fine discretization is used. For the particular selection of model parameters, the amplitude of the reflected wave is approximately 3 × 10 −7 m. As it will be shown later, the accuracy of this coupling deteriorates very fast as the discretization becomes coarser and this approach will only be used for comparison purposes.
To better quantify and compare the accuracy of each coupling case, the total energy within each domain is computed. As the pulse propagates through the bar, energy is transferred from the FE domain (E FE ) to the PD domain (E PD ) and finally back to the FE domain. Figure 6 illustrates the energy transferring between the domains using now the second coupling case with a modulated incident pulse as an example. The pulse crosses the FE-PD interface at t = 1.19 × 10 −4 s and the energy is transferred from the FE to the PD domain. At t = 2.03 × 10 −4 s, the pulse crosses the PD-FE interface and the energy is transferred back to the FE domain. For comparison, the total energy in an equivalent FE only model is also plotted with dashed line.
If the coupling is adequately accurate, the spurious reflections generated should be minimal and all of the energy should be transferred back to the FE domain. The accuracy of the coupling method is tested by comparing the amount of energy that was transferred to the right-hand side FE model compared to the initially provided energy. The percent error between the two energies (energy drift) is used in the following as a measure of accuracy. Figure 7 depicts the energy drift when the second coupling case is used for different values of the PD horizon and discretization length x PD for both the modulated and the Gaussian incident pulses. As increases, much finer grids are required to achieve comparable accuracy. It is noted that, for macroscale problems, the horizon is typically set to = 3 x PD , as higher values can induce excessive dispersion and increase the required computational time, whereas lower values will lead to grid dependent crack paths. 5,21,46 As x PD and increases, total reflection of the incident wave is observed (energy drift = 1 in Figure 7). This is attributed to the cut-off frequency of the discrete PD model being well below the frequency content of the pulse. In such cases, the pulse never enters the PD domain and is reflected back to the FE domain. The parameters x PD and are of outmost importance as they control the model's cut-off frequency.
When the third coupling case is implemented, two additional parameters need to be taken into consideration, namely, the overlapping length l coupl and the choice of blending function h(x) implemented within the formulation. In the work of Aubertin et al, 16 numerical tests were carried out using a linear blending function. Here, a cubic blending function is also used. In the works of Guidault and Belytschko 44 and Bauman et al, 45 a more detailed analysis on the selection of blending functions can be found for the Arlequin method. The problem parameters are the same as before and the number of Lagrange multipliers in the coupling region is the same as the corresponding PD particles.
Using the original PD formulation and the modulated pulse, a parametric investigation is performed for different values of x PD , , and l coupl . In Figure 8, waterfall plots of the resulting energy drift are presented. In all cases, the conditions for total wave reflection remain unaffected by the choice of blending function and l coupl . However, higher values of l coupl can improve the accuracy of the coupling. Furthermore, by implementing a cubic blending function instead of a linear in the formulation, the accuracy of the coupling is improved for the same values of x PD , , and l coupl . This is intuitively expected since the transition is smoother in the second case. Increasing the value of l coupl increases the computational cost of the method as additional PD particles and Lagrange multipliers are introduced. For the cases considered, an overlapping length of approximately l coupl = 0.01 m seems a fair compromise between accuracy and computational cost.
Using l coupl = 0.01 m and a cubic blending function, the energy drift is also computed for the Gaussian pulse and the results are plotted in Figure 9. Comparing Figure 7 and Figure 9, both coupling approaches exhibit total wave reflection under the same conditions, indicating that the cut-off frequency is not affected by the coupling strategy adopted. In both cases, the higher frequency content of the modulated pulse requires refinement of the PD domain to ensure transmission. Furthermore, the accuracy achieved is comparable between the two approaches for the specific choice of parameters. If the linear blending was implemented further refinement would be necessary.
To make this investigation more general, it is of interest to continue the analysis in the frequency domain. To this end, in the next section, the linearized PD theory is coupled with a conntinuum instead of an FE approximation. Furthermore, a brief derrivation of the PD dispersion curve from the discretized equation of motion is presented as it provides insights on the cut-off frequency and its dependence on x PD and that is pivotal in minimizing the spurious reflections.

Dispersion curve for bond-based PD
According to the theory of elasticity, the phase velocity of a travelling wave is independent of its angular frequency and wavenumber , leading to a linear relation between and . In discretized systems however, this relation is nonlinear. Along the interface of two models with different dispersion relationships, spurious reflections occur as the incident wave travels with different speeds in each domain. 47 In PD, waves of higher frequency will travel slower than lower frequency waves. This effect needs to be taken into consideration when coupling domains that are described by different models. Additionally, discrete models exhibit a cut-off frequency that can lead to full reflection of the incident pulse. In the works of Silling et al 28,48 and Dayal, 49 the PD dispersion curve was extracted by substituting solutions with the form u = Ue i( x − t) into the linearized PD equation of motion and solving the resulting integral in the eigenvalue problem. In the work of Mikata, 50 the dispersion relation was derived for different micromoduli, whereas, in the work of Diyaroglu et al, 51 they are given for PD beams and plates. The dispersive characteristic of the PD theory can be also used to match the experimentally observed dispersion, as presented in the work of Butt et al. 52 Considering a bar of infinite length and using the idealisation that the bar consists of a repeating series of masses and springs, the dispersion curve of linearized PD is derived after substitution of Equation (10) in the discretized equation of motion (11). If each particle interacts only with those that is in direct contact with = x and x is uniform, then the equation of motion for a particle becomesü where the volume correction is also included. When u = Ue i( t+ x) , the following relationship is extracted: In the limit when x is very small, we can use the approximation sin( x/2) = x/2 and Equation (38) leads to / ≈ c. Thus, the PD dispersion relation approximates the one extracted from the local theory of elasticity. This observation has already been described in the work of Silling. 28 Following the same procedure, the dispersion relation for various horizon lengths can be derived. As a general case, with = n x, the dispersion relation can be written as follows: where G n (x) is a function whose expression depends on the value of the PD horizon. Expressions for G n (x) are given in Table 1, for n = / x = 1,2,3, and 4. As expected, for x → 0, then G n (x) → n and Equation (39) reduces to the dispersion relationship obtained from classical elasticity.  A comparison between the dispersion relation for different using Equation (39) is illustrated in Figure 10. In the same figure, the linear elasticity case has been included. It is noted that, when the value of the PD horizon is = x, the relationship coincides with that for a two-node linear FE. The figure is plotted for wave numbers ∈ (0, / x) representing only the positive half of the first Brillouin zone. The parameters assumed are E = 200 GPa, = 7850 kg/m 3 , x = 2 × 10 −4 m, and l = 0.4 m. An important feature for this study is that the PD cut-off frequency depends not only on the discretization length x but also on the value of the horizon . In Figure 10, the angular frequencies 0 , which define the cut-off frequency, are indicated for each value of . Values have been nondimensionalised with respect to c = √ E∕ and l.

Continuum-PD coupling and frequency domain analyses
Coupling of linearized PD theory with continuum allows the use of frequency domain analyses assuming rightward harmonic pulses propagating in an infinite bar. Domain (see Figure 11) is divided into two subdomains PD and cont , with PD ∪ cont = and PD ∩ cont = coupl . The classical theory of elasticity will be used in to describe the behaviour of the cont domain, whereas the linearized bond-based PD theory will be used for PD . Depending on the coupling approach, the coupling domain coupl can contain a single material position or an overlapping length. The PD theory is used in the middle of the bar with PD = [−l, l], whereas the rest of the bar is modelled using continuum elasticity with cont = cont,1 ∪ cont,2 = (−∞, −l] ∪ [l, +∞). The incident pulse propagates from −∞ toward ∞. Hence, in cont,1 , there are two pulses, ie, the initial incident pulse and the reflection from the first interface. In cont,2 , on the other hand, only the finally transmitted pulse exists. Denoting with u R , u T , and d the displacement fields in cont,1 , cont,2 , and PD , respectively, we assume solutions in the form u R = U R e −i t , u T = U T e −i t , and d = De −i t , and the governing equations become where U R and U T are solutions in the form U R = e + Re − and (41a) with R and T being the reflection and transmission coefficients, respectively. 53 Energy conservation implies that E R + E T = 1, where E R = |R| 2 is the reflected energy and E T = |T| 2 is the transmitted energy. Using these definitions, the ability of a pulse to be transmitted across a coupling interface is evaluated. Realisation of the first and second coupling case for the coupling of the discretized PD equation of motion with the continuum solution is carried out following the procedure presented in Sections 3.1 and 3.2, respectively. The first case requires the position of PD particles at the interface positions x 1 = −l and x n = l, where n is the total number of PD particles. Continuity and force equilibrium requirements are enforced at the interface positions resulting in the following equations: Using as an example the case where = 2 x, Equations (42c) and (42d) become ( where, earlier, D 1 , D 2 , D 3 , D n − 2 , D n − 1 , and D n are the particle displacements corresponding to the particles with positions x = − l, − l + x, − l + 2 x, l − 2 x, l − x, and l, respectively. Combining Equations (40a) to (40c), (41a) and (41b), (42a) to (42d), and (43a) and (43b), the system of equations can be solved for the R and T coefficients.
In the same manner, the system of equations for the second coupling case can be constructed. Without loss of generality, the example = 2 x is implemented leading to the introduction of two ghost particles at each side of the PD domain. The continuity and force equilibrium requirements lead now to the following equations: where x g1 , x g2 , x g3 and x g4 are the positions of the ghost particles corresponding to the locations x = − l − x/2, − l − 3 x/2, l + x/2, and l+3 x/2, respectively. Although the force equilibrium of the second coupling case, described in Equations (44e) to (44f), looks similar to that of the first coupling case, described in Equations (42c) to (42d), it correlates the sum of the internal forces of each bond that crosses the interface with the internal forces from elasticity. For the example considered here, Equations (44e) to (44f) are expressed as follows:  For the third coupling case to be implemented, PD is extended within cont,1 and cont,2 by a length l coupl = x B − x A = x D − x C to create the overlapping regions coupl,1 and coupl,2 , over which the coupling scheme is enforced ( Figure 12). The assumption of solutions in the form given by Equations (41a) and (41b) is no longer valid within these domains. Assuming a linear blending function a(x) = 1 − (x − x A )/l coupl , an analytical derivation of the unknown fields U 1 (x) and U 2 (x) is sought. Using Equation (30a) and substituting u 1 = e −i t U 1 (x) in coupl,1 , the following differential equation is written: where (·) is the Dirac delta function. A solution in the form j = e i t j (x) was assumed for the Lagrange multipliers and n and x are the total number and the positions of the Lagrange multipliers. A similar equation can be written for coupl,2 .
Using the change of variables z = a(x), Equation (46) becomes where b 1 = √ 2 l 2 ∕E and b 2 = l/EA are two constants andx = ( ∕l coupl . This differential equation can be solved using the Hankel transform. 54,55 The solution obtained in coupl,1 is given as follows: where C 1 and C 2 are two constants, J k and Y k are the Bessel functions of the first and second kind with order k, and H(·) is the Heaviside function. As can be seen from Equation (48), as z → 0 (ie, x → x B ), This solution is not accepted as a finite result is expected and hence C 2 = 0. On the interface between coupl,1 and cont,1 , continuity and equilibrium conditions need to be satisfied. This leads to the following two equations: and Similarly, the same procedure can be repeated in coupl,2 to relate the transmission coefficient T with the solution U 2 . Using Equations (30b) and (30c), (49), and (50), the system of equations can be solved to compute the reflected and transmitted energy across the different domains. The accuracy of the third coupling approach can be improved using a cubic (Hermitian) blending function with a (x) = 3 in Equation (46). In this case, the unknown displacement fields are approximated numerically

FIGURE 15
Reflection -Transmission diagrams for coupling cases 1, 2, 3 with a linear blending function (LB) and 3 with a cubic blending function (CB). The PD discretization was set to x PD = 2 × 10 −4 m. Each row corresponds to a PD horizon value given as = n x PD All coupling approaches considered can achieve accurate wave propagation with minimal reflections for relatively small values of , as |R| 2 is close to 0 and |T| 2 = 1. As the frequency increases, the accuracy of the coupling deteriorates. When the forcing frequency of the pulse is close to or exceeds the cut-off frequency total reflection occurs and the transmitted energy vanishes. The PD discretization length x PD implemented in Figure 15 is the same as that in Figure 10. For all cases, total reflection occurs at the frequency value predicted from the PD dispersion curve regardless of the coupling method chosen. Furthermore, the influence of x PD and on the model's cut-off frequency is illustrated.
It is evident that, as the frequency of the input pulse approaches the model's cut-off frequency, the accuracy of the first coupling case deteriorates faster compared to the other methods. This means that further refinement is required to minimise the spurious reflections at the interface. The values of the reflected |R| 2 and transmitted energy |T| 2 oscillate till the cut-off frequency where total reflection occurs. Close-ups illustrating this fluctuation are provided in both Figure 14 and Figure 15. In this zone, the transmission coefficient is high for certain frequencies indicating that, if the problem solved contains monochromatic incident pulses, the method could be tailored to allow for accurate transmission. However, if the incident pulse excites a broad range of frequencies, only the frequencies corresponding to peaks will transmit. Using the other coupling cases, this zone is much smaller and restricted very close to the cut-off frequency, allowing to accurately simulate wave propagation with coarser grids for the same frequencies.
Using the results reported in Figure 14 and Figure 15, the third coupling case with the cubic blending function leads to the best results, followed by the second coupling case. Through these numerical tests, the improvement of implementing a cubic blending function within the formulation of the third coupling case is also highlighted as the coupling remains accurate for values closer to the cut-off frequency. It is noted here that, for the results presented, the overlapping length was set to l coupl = 10 x. The coupling could be further improved by increasing the overlapping length, as illustrated in Figure 16. The accuracy of the coupling is improved, even for frequencies very close to the cut-off frequency with additional computational cost though.
Finally, of major importance is the implementation of the surface correction factor (Figure 17). The third coupling case is less affected from the skin effect, which appears at the particles near the boundary of PD , where the value of the blending function is almost zero (ie, energy has been transferred to cont ). The first two coupling approaches are affected though, as the coupling interface is located near the area where the skin effect manifests. Comparing Figure 17 with the  Figure 15, the accuracy of the coupling is affected significantly if appropriate correction is not implemented. The third coupling case can lead to the most accurate coupling between PD and continua, given that an appropriate blending function and an adequately large overlapping region are provided. This increases the computational burden of the method and creates ambiguities on the a priori selection of blending functions in 2D and 3D problems. On the other hand, coupling approaches similar to the second coupling case or the method presented in the work of Zaccariotto et al 7 achieve comparable accuracy with lower cost and are easier to implement. The simplistic coupling presented in the first method should be avoided not only because it leads to strong reflections for relatively small frequencies but also because it requires conforming meshes that are problematic in problems involving higher dimensions.

EXTENSION OF THE SECOND COUPLING CASE TO 2D PROBLEMS
The second coupling case is extended to 2D problems to further illustrate the importance of considering the cut-off frequencies of both the PD and FE models when addressing dynamic problems. The only thing that differs in the formulation of the second coupling case is the computation of the forces applied on the interface. In the 1D bar example, the computation of PD interfacial forces required to simply identify which bonds cross the interface (ie, bonds that connect a ghost and a normal particle) and what is the total force of these bonds. In 2D, the interface between the two descriptions is no longer a point but a line and thus it is necessary to compute the exact location where each bond crosses the interface. This is the location where PD forces are applied on the interface. In the following examples, linearized PD are coupled with four-node linear elasticity FE. Assume that nodes 1 and 4, illustrated in Figure 18, define the PD-FE interface and particle j is a ghost particle. From the definition of linearized PD, the bond forces, with respect to the global axes, between particles i and j are computed as follows: where x d j,y } T is the vector containing the displacements of particles i and j, R is the transformation matrix, and K is the local stiffness matrix in 2D given as follows: The value k bond i, is computed using Equation (22). From Equation (8), it is implied that f i,j is equal and collinear to f j,i . In 2D, it is convenient to treat the PD bond forces as externally applied forces to the FE that define the interface. Contrary to that, in the 1D case, the stiffness matrix of the coupled model was constructed directly in an explicit form by adding stiffness contributions. The force f inter m = − f i, = f ,i is applied on the FE and it is concentrated at the intersection point (x A , y A ) between the bond and the coupling interface and subscript m denotes the bond number. The point of application is defined using two Dirac delta functions. Assuming there are no other externally applied forces at the FE domain, then where N are the shape functions of the four-node element and m tot is the total number of bonds that cross the coupling interface. Since the shape functions for nodes 2,3 vanish at the force application location, the interface force is effectively distributed between the two nodes that lie on the coupling interface. The final vector of external forces is computed by summing all the interface forces applied on the interface. Using matrix notation, the equation of motion in the FE domain can be written as follows: The displacements of the ghost particles are again computed by interpolating the nodal displacements of their element using its shape functions as follows: The equation of motion for the particles in the PD domain is Equation (11). Now, the FE and PD equations are solved separately but are coupled together through Equation (55) and the external forcing in Equation (54). Time marching is performed again using the central difference method defined in Equations (16a) and (16b) and (17a) and (17b).
In the 1D problems, the bond constant k( ) was assumed to be uniform within the PD horizon. For the 2D cases, a conical variation of the bond constant is implemented as it provides better convergence to classical elasticity. 40 Assuming plane stress conditions, the bond constant is expressed as follows:

Verification of the 2D model
To validate the coupling in 2D, a pulse propagating in a solid under plane stress conditions is considered. The boundary conditions and geometry are illustrated in Figure 19   The coupled problem was discretized with x FE = y FE = 2 × 10 −3 m, and x PD = y PD = 2.07 × 10 −4 m. The length of the PD horizon was set to = 5 x PD . In total, 123 492 particles were used in the PD domain, out of which, 7330 where ghost particles. The FE element domain was discretized using 30 000 elements corresponding to 30 450 nodes. The total simulation time is t tot = 1 × 10 −4 s with steps t = 5 × 10 −8 s. The standard FE solution was obtained using a uniform grid with x FE = y FE = 2 × m. In Figure 19 (right), an example of the discretization of FE and PD is illustrated. The location of the PD particles is independent of the FE nodal locations, thus simplifying the process. If the first coupling case was to be implemented, refinement of the FE mesh would be required near the interface for the FE nodes to coincide with the PD particles.
In the left column of Figure 20, the propagation of the Gaussian pulse is illustrated at time instants t = 4.4 × 10 −5 s, t = 5.7 × 10 −5 s, and t = 9.9 × 10 −5 s. Similar to the 1D case, each time the pulse crosses a PD-FE interface, a reflection is generated. Since a much finer discretization was implemented in the PD domain compared to the FE domain, the cut-off frequency of the PD domain was higher than that of the FE domain, allowing the pulse to be transferred. This can also be seen in the middle column of Figure 20, where the displacements along the reference line, illustrated in Figure 19, are compared for the coupled and the FE-only solutions. Compared to the amplitude of the pulse, the amplitude of spurious reflections is negligible, and the two solutions agree well. To further compare the two solutions, the L 2 norm of error is employed. Since the nodal points of the two solutions do not coincide in the PD region, the solution in the PD region is interpolated onto the FE nodal points through linear interpolation. The magnitude of error is defined as follows: whereû coupled contains the interpolated values of the coupled model u coupled onto the FE nodes. The magnitude of the L 2 norm is plotted in the right column of Figure 20. As expected, discrepancies first appear when the pulse crosses for the first time the FE-PD interface and reach the maximum value of 2.02 × 10 −6 when the pulse exits the PD domain. It is noted that the error here is not only due to the generated reflections but also due to differences in the numerical dispersion characteristics of the two solution methods. Still, the two solutions are in very good agreement in the whole computational domain.

Partial trapping of pulses in PD-FE models
According to the findings of the 1D study, if the cut-off frequency of the PD domain in the previous example was not high enough to allow the pulse to be transmitted accurately, the spurious reflections generated would be more severe and affect the accuracy of the simulation. When coupling methods such as PDs with FE, it is desirable that the PD domain will be localised only in a small area of the computational domain, where steep stress gradients are expected. To accurately capture these phenomena, very fine discretization is usually required. It is thus expected that, in practical applications, the PD grid will be much finer than the FE mesh. Since the frequency content of the pulse was low enough to be transmitted within the coarse FE description, it was also able to propagate within the PD domain. In this case, consideration of the dispersion curves for each model is needed to make sure that the cut-off frequency of the FE domain is high enough to allow the pulse to be transmitted. At the same time, the cut-off frequency of the PD domain needs to be higher to reduce the spurious reflections due to the coupling. Consideration of the FE cut-off frequency is of equal importance to the PD one to achieve accurate coupling. This is evident in cases where the source of the pulse is located within the PD domain. One such case is considered in the next example.
Assuming an infinite plate with a hole at its center, the area around the hole is modelled with PD, whereas the rest of the plate with FE. To avoid interference due to the boundary conditions of the numerical model, the computational domain is truncated by adding a perfectly matched layer (PML) to attenuate outward travelling waves. In a recent study, unbounded problems were addressed using the PD differential operator and Sommerfeld boundary conditions. 56 The PML was originally proposed by Berenger for applications in electromagnetic waves 57 and was later extended to applications involving elastic waves. 58 Here, the convolution PML (C-PML) proposed in the work of Matzen 59 for time-domain analyses is employed. In the PML region, the coordinate variables x i , i = 1,2 of the problem are transformed into the stretched coordinate variablesx i , in each direction, given as 58,60 where is the angular frequency and s i are the stretched coordinate functions, proposed by Kuzuoglu and Mittra, 61 ie, Three coordinate-wise functions are introduced in Equation (59), i (x i ) ≥ 0, (x i , ) ≥ 1, and i (x i ) ≥ 0. i (x i ) controls the attenuation in the PML region, whereas (x i , ) and i (x i ) enhance the attenuation of evanescent waves. In the work of Matzen, 59 it was shown that, when near-grazing waves do not manifest, the selection (x i , ) = 1 and i (x i ) = 0 simplifies the PML computations as one of the convolution terms disappears without affecting significantly the accuracy. Since this is applicable to the case considered here, these values are adopted. The function i (x i ) is defined in the works of Matzen 59 and Madsen et al 60 and the references therein as follows: where max i = − ( 4c p log 10 R 0 ) ∕ (2d i ), R 0 = 10 −8 is the theoretical reflection coefficient at normal incidence, x 0 i is the location of the PML interface, and d i is the PML thickness in each direction. The coordinate stretching performed in the PML region leads to a nonhomogeneous differential operator. To extract the original one, artificial anisotropy is introduced to the material. The FE equation of motion now is where M, Z, K are the mass, damping, and stiffness matrices, and g is a convolution term defined in. 53 The interested reader can find more details in the work of Ben-Menahem and Singh 53 on the FE implementation of the PML, whereas a PML formulation for state-based PDs can be found in the work of Wildman and Gazonas. 62 The geometry of the problem is illustrated in Figure 21 along with the definition of the coordinate-wise function i in different areas of the plate. A radial pressure is applied at the walls of the hole with F (t) = F max e . The problem solution is repeated two times; the first time the FE discretization was refined to minimise the spurious reflections generated, whereas, for the second time, a coarser FE mesh was used. In both solutions, the PD and the FE discretizations in the PML were not changed. The model parameters are summarized in Table 2. The displacements are constrained on the outer perimeter of the PML as according to the work of Matzen, 59 this improves the stability during time marching.
The results using the coarser and finer discretization are illustrated in Figure 22 and Figure 23. The displacement magnitude |U| = √ U 2 x + U 2 is plotted at different time instants as the pulse propagates along the PD-FE and FE-PML interfaces till it is finally attenuated. In the same figures, the total energy in the PD and the FE model is also included. When the fine FE discretization is used, the pulse propagates without significant reflections across the PD-FE interface. The total energy provided into the PD model is transferred to the FE model and subsequently exits the domain (vanishes due to the introduction of the PML). On the other hand, when the coarse mesh is used, only part of the energy is transmitted to the FE domain, whereas approximately 16% of the energy is trapped within the PD domain.
By making the FE mesh coarser, the cut-off frequency is reduced, and significant spurious reflections are generated at the interface. These reflections are trapped in the PD domain and are unable to reach the PML region. In practice,

FIGURE 22
Pulse propagation with close up in PD and energy transferring using the fine FE mesh

FIGURE 23
Pulse propagation with close up in PD and energy transferring using the coarse FE mesh. Approximately 16% of the total energy is trapped within the peridynamic domain if the frequency content of the propagating pulse is known beforehand, the discretization parameters of the numerical descriptions will be selected accordingly. However, application of impact-like loads in the PD domain will excite the whole frequency spectrum of the finer model. It is thus desired that the two models have similar cut-off frequencies to minimise the reflected energy. In problems involving dynamic crack propagation, mismatch between the two descriptions will lead to accumulation of the trapped energy within the PD domain that will lead to erroneous and unrealistic results.

CONCLUSIONS
Three approaches to couple FE meshes with PD grids have been formulated to study the spurious reflections that appear during dynamic simulations. Although these reflections cannot be avoided completely, they can be minimised through careful selection of the model parameters. The first coupling approach presented is the simplest one and used mainly for comparison purposes. Its use is restrictive as it requires conforming meshes, posing difficulties for 2D and 3D applications. 14 The second and third approach can be extended to problems with higher dimensions more efficiently, 7,26,32,63 as both of them can couple models with no restrictions on the position of the PD particles and the FE nodes.
The coupling approaches were compared through numerical tests in a 1D bar for pulses with different shapes and frequency contents. Initially, the energy drift is computed in the time domain for incident pulses with modular and Gaussian shape. The numerical tests are then continued in the frequency domain where the discrete PD model is coupled with continuum mechanics. This allows to study more systematically the accuracy of each coupling approach by computing the transmission and the reflection coefficients in an infinite bar. The discretization size and the value of the PD horizon are of outmost importance toward an efficient coupling as they control the cut-off frequency of the numerical solution.
The results indicate that the third coupling case with a cubic blending function leads to the most accurate energy transmission. Comparable accuracy can be achieved implementing the second coupling case (or similar methods presented in the literature, eg, the work of Zaccariotto et al 7 ). This approach is simpler to formulate, less computationally expensive, and does not entail ambiguities with regards to the optimal selection and construction of blending functions.
The second coupling case is also extended to 2D problems to illustrate the significance of considering the cut-off frequencies of both the finer and the coarser models. An axially loaded hole in an infinite domain is used as an example. Analytical solutions of similar problems can be found in the work of Graff. 64 In dynamic problems, where the source of the pulse is located within the fine PD model (eg, dynamic crack propagation), entrapment of the energy can be observed. In such situations, cut-off frequency of the coarser model needs to be adjusted to match that of the finer model.