Imposing local boundary conditions in nonlocal/peridynamic models is often desired/needed. Fictitious nodes methods (FNMs) are commonly used techniques for this purpose but they are limited, in general, to domains with simple geometry. FNMs also mitigate the well-known peridynamic surface/skin effect at boundaries/surfaces. Here, we introduce a general algorithm that automatically locates mirror nodes for fictitious nodes in the mirror-based FNM, without requiring an explicit mathematical description of the boundary. The algorithm is based on computing a nonlocal gradient, at fictitious nodes, to determine the “generalized” normal direction to the boundary of a domain with arbitrary geometry. We test several FNMs on peridynamic diffusion problems with or without singularities, that exist in the corresponding local models, along the boundaries. We find that the mirror-based FNM works best in agreeing with the classical solutions. We then test the new algorithm with this FNM for diffusion problems in domains with complex geometries, including one with intersecting cracks. The algorithm is general and should work for any type of peridynamic model, including those for problems with moving boundaries and growing cracks, for which enforcement of local-type boundary conditions is desired. Since high accuracy is critical near boundaries of arbitrary shape (including corners, notches, crack tips) in a variety of problems, the new algorithm has potential for high impact. Keywords Boundary conditions; Diffusion; Fictitious nodes method; Nonlocal gradient; Peridynamics; Surface effect. 1. Introduction Two kinds of singularities usually appear in classical PDEs: angular and interface singularities [1]. Angular singularities occur at reentrant corners (including crack tips), while interface singularity happens due to sudden jumps in a boundary condition (BC). While angular singularities are more commonly observed in models of physical systems, there are certain problems in which both may exist. For example, in galvanic corrosion [2], there is a sudden change of boundary condition at the interface between anode and cathode, and at the same location, a reentrant corner forms due to localized corrosion. They lead to singularities in the classical formulations of the Laplace/Poisson equation which governs the electro-potential distribution * Corresponding author. E-mail address: fbobaru2@unl.edu (F. Bobaru). 2 in the solution domain. Classical methods deal with singularities by special treatments such as conforming transformation method [3], local refinement [4,5], singular function methods [6–8], combined methods [9,10], etc. These methods, while working fine for simple cases, have limitations on problems with moving singularities such as crack propagation and corrosion. The peridynamic (PD) model was introduced by Silling in 2000 [11] and it allows for a natural treatment of discontinuities/singularities by employing integration, over a nonlocal region called the horizon region, rather than differentiation. For a nonlocal formulation, associated BCs are of nonlocal type as well, and they are sometimes referred to as “nonlocal volume constraints” [12,13]. In reality, however, conditions to be imposed (on values and/or derivatives of the unknown function) are known (measurable) only at the surface of a body, not through a finite thickness layer at the surface. The natural representation of such measurement-based conditions is via local boundary conditions. Therefore, imposing local BCs in nonlocal/peridynamic models is often desired/needed. Fictitious nodes methods (FNMs) [14–17] are commonly used to enforce local boundary conditions in PD models. Other names for fictitious nodes are: “dummy particles”, “ghost particles”, or “imaginary particles” [18,19]. When applied to free boundaries (or crack surfaces) they can also mitigate the PD surface effect [20]. The PD surface effect appears because, unlike in the bulk, points near the surface/boundary do not have a full nonlocal neighborhood. One of the FNMs, which uses mirror nodes and thus is called the mirror-based FNM, has been proven to work efficiently in problems with simple geometries [16,17,21]. For arbitrary geometries such as those with curved boundaries, kinks, corners and cracks, however, the method is difficult to implement because there is no general algorithm to find the mirror nodes (of all fictitious nodes) required by the method. In this work, we introduce a new algorithm that helps us automatically find mirror nodes in the mirror- based FNM for domains of arbitrary geometry, including for those with cracks. This algorithm approximates the nonlocal normal vector along the boundary for any geometry by computing the PD gradient of the “index” function. The index function, for a generic point , is the proportion of the area of ’s horizon that is inside the domain and connected to , to the full area (or volume in 3D) of ’s horizon. Note that parts of the horizon region may fall outside of the domain but they are counted in the “full area” of the horizon region. We select the PD diffusion model to test our algorithm. With small modifications, the algorithm presented here is also applicable to PD models solving other types of problems, such as fracture and corrosion damage. We also investigate two other types of FNMs and compare their results with those from the mirror-based FNM: the “naive” FNM [22], and the Taylor FNM [18,23,24]. We compare the performance of different FNMs in enforcing local boundary conditions using two problems: one is a simple problem without singularity and the other is the Motz’s problem with a singularity (in fluxes) along a boundary where Dirichlet and Neumann boundary conditions meet [3,25]. We show the capability of the new algorithm, used in the context of mirror-based FNM, by solving examples of diffusion in domains with crack surfaces and curved boundaries. This paper is organized as follows: in Section 2 we review the PD method for diffusion-type problems; in Section 3 we discuss the fictitious nodes method and introduce the autonomous algorithm to generalize the mirror-based FNM for arbitrary geometries; in Section 4 we compare the performance of three different types of FNMs using examples with and without local singularities, then test the generality and capability of the developed algorithm on the mirror-based FNM for more complicated domains (with cracks); conclusions are given in Section 5. 3 2. The peridynamic model fordiffusion The peridynamic (PD) theory is a nonlocal extension of the classical continuum mechanics [11]. The PD formulation is in the form of integro-differential equations (IDEs) rather than PDEs used in classical local theories. In PD, discontinuities such as cracks/damages can initiate and propagate naturally and autonomously [15,26–28]. While PD has been primarily used to deal with mechanical behaviors [26,27,29– 32], it has also been employed in diffusion-type problems involving cracks and damages, including thermal diffusion [14,33–35] and mass transport (e.g. corrosion) [36–42]. In peridynamics, each spatial point interacts with other points within its neighborhood which is called the horizon region of and is usually selected to be a disk in 2D (or sphere in 3D) centered at . The radius of is called horizon size (or simply, the horizon) and denoted by . Objects that carry the pairwise nonlocal interactions between points are called bonds. Fig. 1 schematically shows a peridynamic body with a generic point , its family and its horizon. Fig. 1. A peridynamic body with a generic point and its horizon . Nonlocal interactions exist through the bond between two points, e.g., point and an arbitrary point located in its horizon . In this work we use the bond-based peridynamics for diffusion-type problems. In the context of heat conduction, for example, one assumes that the heat diffusion between two points depends only on their temperature difference. In the more generalized state-based formulation, heat diffusion in a bond is directly influenced by other bonds which connect to the same point [14]. This setting may be beneficial for certain problems but is not considered here. For diffusion in a homogeneous and isotropic medium with constant diffusivity, the peridynamic transient equation can be written as [12,33,43]: (,)=(,)+(,) (1) where (,) is the target quantity (e.g. temperature, concentration, electric potential, etc.) at point at time , is the scalar diffusivity constant, is the source/sink term and is the PD Laplacian operator which can be expressed as: (,)=()[(,)(,)]d, (2) 4 where the kernel function is nonnegative and radially symmetric, i.e. ()=()0 and d is the differential volume (area in 2D, length in 1D) of material point . It has been shown that as the horizon approaches zero, in a continuous domain, the peridynamic solution to a diffusion problem converges to the solution given by the classical Fick’s law of diffusion, under certain conditions [22,44,45]. A kernel function that leads to good convergence properties with the one-point Gaussian quadrature discretization of the spatial domain is [22]: ()= (3) where is a constant calibrated by matching PD flux to classical flux for a linearly distributed field. Other options for the kernel function (see [22]) can be selected, but they are not considered here. The steady state of Eq. (1) without sources/sinks is: (,)(,)d= 0 (4) For a specific problem, nonlocal boundary conditions (sometimes called “volume-constraints” [12,13]) are required to solve Eq. (4). In the next section we discuss local and nonlocal boundary conditions in peridynamics. 3. Fictitious nodes methods to enforce local-type boundary conditions and mitigate the peridynamic surface effect Unlike classical local methods, the boundary conditions in peridynamics are nonlocal. However, when solving practical problems, imposing local-type boundary conditions in nonlocal/peridynamic models is usually desired/needed because, in reality, conditions (on the unknown function values or its flux) are imposed at the surfaces of a body, not through a finite layer near the surface. The natural representation of such conditions (based on measurements) is via local boundary conditions. Various methods to impose local boundary conditions in PD models have been investigated in [13,14,46]. One such method is the fictitious nodes method (FNM) [14,16,46]. In FNM for peridynamics, instead of imposing local BCs at boundary points, certain volume constraints are specified on fictitiou the “collar” outside of the body shown in Fig. 2, called the fictitious region), so that desired local boundary conditions are satisfied or approximately satisfied. The fictitious region is the smallest region that assures full horizon for all domain points. Such volume-constrained peridynamic problems are defined analogous to boundary value problems with PDEs in the local theory according to nonlocal vector calculus [13]. Volume constrained peridynamic steady-state diffusion (or Laplace) problem can be expressed as [13]: (,)=0 ()= 0 (volume constraints ) (5) where is a function that prescribes the constraints of on . Fig. 2 schematically shows the solution domain , its boundary , and the constrained volume or fictitious region, . 5 Fig. 2. Schematic of a peridynamic domain ( ), its boundary ( ), and its constrained volume (or fictitious region, ). When applied to free boundaries (or crack surfaces), the FNM mitigates the peridynamic surface effect [16]. This is a common issue for nonlocal methods as well as local methods with certain nonlocal treatments such as the SPH method [19,23]. The surface effect appears because, unlike in the bulk, material points near the surface/boundary do not have a full non-local neighborhood. The surface effect leads to slightly different behavior of material points near the surface compared with those in the bulk. This could mean a lower diffusivity near the surface for diffusion problems, and a lower stiffness for problems in elasticity (see [20]). These effects are reduced as one decreases the horizon size and would be “exact” when the size of the neighborhood is the same as the physical nonlocal interaction range, which could be atomistically small. In practical modeling, the size of nonlocal interaction is usually set to match observable physical length-scales (see discussion in [47]), and not larger than relevant geometrical features of the domain (notch widths, etc.). In this section, we briefly review three different types of FNMs from literature and then introduce the new algorithm for the mirror-based FNM to make it applicable to arbitrary geometries. Comparisons between their performances are shown in Section 4. 3.1. Naïve FNM A naïve type of FNM is often used in the literature to impose local Dirichlet and homogeneous Neumann (no flux) BCs. This method enforces Dirichlet BCs by assigning the same values to all fictitious points corresponding to a boundary point, while homogeneous Neumann BCs are enforced by simply neglecting all corresponding fictitious points [22]. See Fig. 3 for an illustration of how a Dirichlet BC (,)= for (boundary subjected to Dirichlet BC) is enforced at each fictitious point . An even simpler version of “enforcing” local Dirichlet BCs in the nonlocal model is to dispense entirely with the fictitious points and subject only the surface points to the values of local boundary conditions [33]. Previous work has shown that results by both versions converge to local solutions as the horizon size approaches zero [22,35]. The naïve FNM has the advantage of featuring the simplest implementation and is the most efficient (see Section 4.1). Note, however, that a jump-discontinuity in the solution at the boundary may be generated, leading to possible errors in fluxes near the boundary (see the result in section 4.1 for an example). 6 Fig. 3. Illustration of enforcing the local Dirichlet BC in the naïveFNM (redrawn from [19]). 3.2. Taylor-based FNM The second FNM, used in the PD context first in [24], requires a Taylor expansion (to linear terms) for points in the horizon of the boundary. We call it Taylor-based FNM or simply Taylor FNM. To impose a local Dirichlet BC (,)=() for (boundary subjected to Dirichlet BC) ,the constraint value at a generic fictitious point , which is located inside the horizon of an generic point in the domain with intersecting the boundary at , is extrapolated from as: =+ (+)+()+( )(+) =1+ (6) where represents (,), represents (,), and represent the distance of point and to the boundary, respectively. Notice that with this approach, the imposed constraint at a fictitious point has to be computed for each corresponding domain point anew, and this means that such value may change during computations at the same solution step (see Fig. 4, where is the distribution of () for , and is different for each ). The solution step refers to each call to the Conjugate Gradient (CG) solver, see Fig. A.2 in Appendix A for the flowchart of the simulation. 7 Fig. 4. Illustration of enforcing a local Dirichlet BC (,)= for with the Taylor FNM (redrawn from [19]). Eq. (6) may lead to unstable results if implemented directly. A modification is required as follows [23]: =+( )×min,1+ (7) where is a parameter determined from numerical tests. Normally =1.5 leads to good results [23]. To impose a local Neumann or Robin BC (,)=(,) for (boundary subjected to Neumann or Robin BC), where is a given function, we can first relate values at a generic domain point and a generic fictitious points (inside the horizon of ) to the BC by the following Taylor expansions to linear terms around (intersection of and ) [24]: =+ +()+() = +() () (8) Note that, since , =(,)is not given and needs to be computed. By subtracting these two equations we have: (+)() (9) where is approximated from the first equation in Eq. (8) and then substituted into Eq. (9) to solve for . Notice that if () is a nonlinear function, a nonlinear equation solver, such as Newton’s method, is needed to solve for [48]. Notice that, if =in Eqs. (6) and (9), the Taylor-based FNM degenerates to the mirror-based FNM which will be introduced in the next section. 8 3.3. Mirror-based FNM Another FNM,referred to as the mirror-basedFNMor simply “mirror FNM”, assignsthe constraint ateach fictitious point based on the current solution of its mirror point in the domain. Two points are a pair of mirror points if they are “symmetric” to the boundary. In this method, one assumes that for each fictitious point , its mirror point ()is given [14] as shown in Fig. 5where () is the mid-point of () sitting on the boundary. When the boundary is smooth in a neighborhood (of size that depends on geometry) of (),()is collinear with the normal vector to the boundary at(). Fig. 5. Mirror-based FNM for imposing local BCs on by assigning certain constraints at fictitious points in based on current values at their mirror points inside , and the value of the boundary condition to be imposed. To enforce the local Dirichlet BC (,)=()for and the Neumann BC (,)=() for in the mirror FNM, the constraint at a generic fictitious point is written as: =2() ()(10) and: =() ()(()) (11) respectively, where =(,), ()=(),, ()=(),,and ()= ()+()()= 2()= 2()() is the distance between and (). The local Robin BC (,)=((,)) for is enforced on fictitious points as: =() ()(())=() () +() 2 (12) in which the approximation ()=() is made according to Eq. (10), and requires to be solved using a nonlinear solver if function is nonlinear. In both the Taylor and the mirror FNMs, constraints applied at fictitious points vary with the solution step (or time step, see the flowchart in Fig. A.2.). However, in the same solution step, the constraint at each fictitious point in the Taylor FNM may change with the domain point at which a computation is performed. 9 For the mirror FNM, this does not change because it only depends on the current solution at its mirror point in the domain (uniquely determined by the algorithm) and the value of the imposed boundary condition (at the intersection point between the boundary and the segment connecting the fictitious point and its mirror point, which is the mid-point on the segment). See Fig. 4 and Fig. 6 for illustrations of how local Dirichlet BCs are enforced in the Taylor and the mirror FNMs, respectively. Fig. 6. Illustration of enforcing a local Dirichlet BC in the mirrorFNM (redrawn from [14]). Important Remark: note that in the original Taylor and mirror FNMs in PD, mirror pairs are considered given as inputs, which works for simple geometries. However, for the general case, to find mirror pairs one needs to find the distance between each point and the boundary, and, for each fictitious point, find the mirror point in the domain. While providing such information is simple and straightforward for simple geometries, it is not so for arbitrary ones, for example shapes with corners of various angles, cusps and crack tips, etc. At such locations, the tangent line (and normal vector) may not be well defined. While some ad-hoc choices can be made for overcoming this issue, we aim for a general algorithm that generates a well-defined “generalized” normal vector to the boundary at any of its points, regular or not, so that both the mirror point and the distance from the point to the boundary can be determined without ambiguity. In the next section we introduce such an automatic algorithm so that local BCs can be enforced in PD models using the mirror FNM for arbitrary geometries. 3.4. An algorithm to find mirror points using the peridynamic gradient Using the mirror FNM for arbitrarily shaped boundaries (with corners, cusps, etc.) requires finding mirror points for all fictitious points. For surfaces of thin cracks, domain points (within a distance from the crack surface) on one side of the crack surface can serve as fictitious points for those on the other side. In what follows, includes all points in the fictitious region as well as domain points (acting as fictitious points) that are required to enforce local boundary conditions. 10 As shown in Fig. 5, ()can be determinedby()=+2(())once(),the mid-point of the segment connecting ()and ,is known. Therefore, we first need to locate ().Starting from , one could determine ()by searching the family points of for the domain point closest to . However, such a domain point may not be unique. In addition, for fictitious points near corners with an angle smaller than 90° (see regions in Fig. 7), this strategy may find mirror points which are not in the domain. Special treatments would have to be used for these instances, which reduces the generality of the algorithm. Fig. 7. Schematic of fictitious regions near an acutecorner for which mirror points to some fictitious points are not located inside the domain, if ()is determined to be the domain point closest to . We introduce a more general approach to compute ()by computing a modified nonlocal gradient (unit) vector at a fictitious point , =(,), which will allow us to locate a unique (). For smooth boundaries, has the same direction with the normal to the boundary at the corresponding location on the boundary. For non-smooth boundaries in the horizon of ()(e.g. at corners, kinks, etc.), may or may not be collinear with the outward normal vector (which may or may not exist) to the domain boundary. To compute in the peridynamic framework, without a mathematical representation of the boundary, we first compute the peridynamic gradient (or flux) [34], at a fictitious point, defined by: (,)= [(,)(,)]d (13) where is ’s horizon, and (,) is the index function given by: (,)=(,)d d(14) where (,) is the binary connectivity function which equals 1 if is connected to , and 0 otherwise. Thus, (,) is the ratio of the area of that is inside the domain and connected to , to the total area (or volume in 3D) of . Although (,) is dependent on time, it changes only when the domain changes. See Fig. 8 for an illustration of (,). 11 Fig. 8. Index function (,)for points near the boundary. Notice that the function does not have to belinear as shown in the picture. We note that (,) as crosses a thin crack is an inverted hat function, instead of the monotonous function shown in Fig. 8 for a generic boundary. In such a case, its peridynamic gradient formulated by Eq. (13) will not guarantee a direction equivalent to the outward normal to the crack surface, which is what we want to have it represent. Modifications of Eq. (13) are then required to resolve this issue as follows: (,)= (,)d (15) where = . For domain points (acting as fictitious points) within a distance from a crack surface, is the part of that is disconnected from .Then at a fictitious point is obtained by: =(,) (,)(16) Notice that has the opposite sign compared with the external normal at a regular point on the boundary. With computed, ()can be determined as the domain point closest to along and ()by +2(()), as shown in Fig. 9. Fig. 9. Schematic diagram for a fictitious point () and its corresponding projection (()) onto the boundary and mirror point(()) in the domain, along direction. 12 Notice that this algorithm may encounter a problem when finding mirror points at sharp (relative to the horizon size) convex corners, but this only affects some very small portions of the fictitious region (see the green regions in Fig. 10). The issue can be reduced by using smaller horizon sizes but not fully resolved. Approximation of the constraint imposed at such points can be achieved by averaging the constraint of surrounding fictitious points. Another option is to use Naïve FNM for such nodes. A comparison between Fig. 7 and Fig. 10 shows the improvement in this situation when we use the more general normal direction (based on PD gradients) instead of the regular normal to the boundary surfaces. Fig. 10. Schematic of fictitious regions near a convex corner for which mirror points cannot be located appropriately, if ()is determined by the domain point closest to along . When implementing the algorithm in the discretized configuration (see Appendix A), as shown in Fig. 11, we first search all the family nodes of for the closest domain node () to (or first node with ,()= 0 for “fictitious” nodes at crack surface) while the distance between ()and is smaller than /2. Notice that if a node sits on the boundary, it belongs to the domain node. In certain cases where two such ()exist, selecting either one of them is acceptable. Then we have ()+ (() /2)and ()=() /2. Ideally, ()==+ 2(), but since may not coincide with a node’s coordinates, instead, we search the domain node closest to . This implementation works for both uniform and non-uniform grids and () converges to the analytical , for a certain choice of , as / approaches infinity. Since the number of fictitious nodes usually only accounts for a small portion of the total number of nodes, the computational cost to locate mirror nodes for all fictitious nodes is trivial compared with the cost of a complete simulation. 13 Fig. 11. Skematic diagram for determing () for a generic with four different possible . The summary of the algorithm’s implementation is: At time , and after any time step at which geometry changes, do: 1) For all nodes in the domain: compute the index function (,) according to Eq. (14). 2) For all fictitious nodes : a. calculate the “generalized” gradient according to Eq. (15). b. compute distances between its family nodes in the domain, and direction; select only those with this distance smaller than /2; from this set of nodes, choose the node (), closest to . c. compute =+(2() ), an approximate point to the exact mirror node for ; find (), the node grid in the domain closest to , by searching in a box near it. The autonomous algorithm developed works for arbitrary geometries, including corners with different angles, cusps and crack tips, where the peridynamic gradient (or ) may not be collinear with . In the original mirror FNM, one treats fictitious nodes at sharp corners either by assuming to be the average of the normal vector on its two edges [14], or assuming to be the same as on one edge [16]. By using the new algorithm, however, the mirror nodes for fictitious nodes at the corner are found automatically and consistently just like the corner is a special case of a curved boundary. As shown in Fig. 12 where those arrows start from fictitious nodes and end at their mirror nodes, has smooth transitions between the normal vector on two edges of the corner (or crack tip). Not all the arrows are shown in order to illustrate clearly how transits. The transition zone will narrow down as the horizon shrinks ( -convergence) and 14 the transition will be smoother as the ratio of and increases ( m-convergence). It may be possible that this transition of the peridynamic gradient be used to detect/track crack tips, on the fly, but this idea will be pursued elsewhere. (a) (b) (c) Fig. 12. Ficititious nodes and their mirror nodes(connected by arrows), computed by the new algorithm, for three diffferent configurations: (a) at a corner; (b) at a crack tip; (c) at the edge of a crack. The algorithm to find mirror nodes can be implemented as a preprocessing subroutine, because mirror nodes do not change if the boundary or crack surface remains unchanged. However, if the solution domain evolves, which happens in problems such as corrosion or fracture, re-finding of mirror nodes is required. Explicit formulas for boundaries or crack surfaces are not required for this algorithm. This algorithm also works for 15 damage without explicit surfaces, such as distributed damage (a subdomain containing nodes with damage indices smaller than the threshold for a crack surface, e.g. 0.4) [49]). This topic, however, is left for future research. For problems with complex shapes, the algorithm leads to important time-savings compared to manually determining mirror nodes. In the following section, the performance of the mirror FNM will be compared with the naïve and the Taylor FNMs. We will also show the capability of the autonomous mirror FNM using problems with curved geometry and cracks. 4. Results and discussion In this section, we first show how previously discussed three different FNMs perform in enforcing local BCs through the steady-state diffusion in a square domain. Then we will show how these FNMs perform in a problem with a singular point on its boundary. Finally, we are going to show the capability of the autonomous mirror FNM in problems with curved geometry and cracks. The classical PDE-based formulation for steady-state diffusion problem can be written as: where (()) defines the boundary conditions (which could be Dirichlet, Neumann, Robin, or mixed). For 2D problems, = (,)in Cartesian coordinates and = (,)in polar coordinates. 4.1. The performance of the three FNM versions in enforcing local BCs in peridynamics for problems without singularities In the first problem, we consider a square domain with its side equal to 0.1, subject to the following local boundary conditions: as shown in Fig. 13. The classical solution to this problem is a linear function written as: In the PD simulations for this problem, the horizon and grid size are taken to be 0.004 and 0.001, respectively. This leads to about 10,000 nodes for a uniform discretization. ()=0 ()= 0 (17) (=0.05,)=0.1 (= 0.05, )= 0.6 (,= ±0.05) = 0 (18) (,)=5+0.35 (19) 16 Fig. 13. The geometry and boundary conditions for diffusion in a square domain used to compare different FNMs’ capability in enforcing local BCs in the PD formulation. The solution, obtained with the different FNMs, along the dashed mid-line in Fig. 13 is shown in Fig. 14. All three results match the classical solution very well. However, if we zoom in near the boundary, we can see clearly that the result obtained with the mirror FNM matches the classical analytical solution much better than results obtained with the other two FNMs. The simulation times for the three FNMs are given in Table 1. The simulation with the naïve FNM is at least 10 times faster than other two types of FNM, while the efficiency of Taylor and the mirror FNMs are similar to one another. Therefore, for such a simple problem, if the accuracy of the solution near the boundary is not critical, the naïve FNM could be taken as the “best” option. If high accuracy is needed, then the mirror FNM should be used. Fig. 14. The classical analytical solution and PD solutions with different FNMs for the steady-state diffusion problem in Fig. 13, along the dashed line shown there. Table 1 Computational time of the simulation using the three different FNMs, when = 0.004 and = 0.001. FNM naïve Taylor-based mirror-based Time (s) 2.0 20.1 20.6 17 Notice that in [24], the Taylor FNM shows slightly better results than the mirror FNM. This is because for the problem studied in [24], Eq. (7) is not required since there are no stability issues in that case. However, for the example shown in Fig. 13, as well as the next example shown below, the solution obtained with the Taylor FNM does not converge if Eq. (7) is not used. This issue of numerical stability depends on the problem being solved. The mirror FNM, however, does not need such extra treatments for numerical stability. The next example has a singular point along a straight boundary, and we will use it to further test the performance of different FNMs types. 4.2. The performance of the three FNMs for Motz’s problem To test the capability of the different FNMs in handling local singularities along the boundary, we choose Motz’s problem [3,25], which can be seen as a steady-state diffusion problem with the following local boundary conditions: as shown in Fig. 15. This classical solution for this problem has a strong local singularity (/) at the origin . Fig. 15. Domain and boundary conditions for an example of Motz’s problem. The classical solution for Motz’s problem can be written as [3]: where ’s are analogous to the stress intensity factors in linear elastic fracture mechanics, sometimes called “generalized flux intensity factors” [50]. Here we choose 34 terms, as in [51]. In Fig. 16 (a) we show the contours for the classical solution. The point-wise relative differences between the PD solutions and the classical solution are given in Fig. 16 (b)-(d). For the PD simulations, we used =0.04 and = 0.01, respectively. Notice that for the mirror FNM, the non-trivial relative difference (>5%) is restricted to the horizon region of the singular point (point ) and the left-bottom corner (point D). As the horizon size approaches zero (with / not decreasing), these areas also converges to zero. For the other two types of FNMs, the relative difference is large not only near the singular point or the corner, but also at locations far away from them. A quantitative comparison of the solutions along the vertical dash (=1,)=0 (,= 1)= 0 (> 0,= 0)= 0 (< 0,= 0)= 0 (= 1,)= 500 (20) (,)=/cos+1 2(21) 18 line shown in Fig. 15 between the three types of FNMs is shown in Fig. 17. As we can see in the zoomed- in image, the PD solution obtained with the mirror FNM matches the classical solution much better than the PD solutions obtained with other two FNMs. (a) (b) (c) (d) Fig. 16. (a) Contours for the classical solution for Motz’s problem; (b)-(d) Relative difference(Rel. Diff.) to classical solution of Motz’s problem using the PD model with the naïve, the Taylor and the mirror FNMs, respectively. The -convergence [52] (with / fixed to be 4) to the classical solution for the three PD solutions at point P (-0.5, -0.48) (see Fig. 15), is provided in Fig. 18 with the relative difference in log scale. We do not choose a point on the boundary because the analytical solution there is zero (i.e. relative difference does not exist). Notice that as the horizon size changes, the grid also changes and there may not exist a node at Fig. 17. The classical solution and PD solutions with different FNMs for Motz’s problem along the vertical dashed line at =/2+/2, as shown in Fig. 15. 19 point P. In such cases, we simply obtain the value at point P by averaging the value at four nearest nodes around P. Fig. 18 demonstrates that, as the horizon size decreases, solutions obtained with all three types of FNM approach the classical solution. However, the mirror FNM produces relative differences from the classical solution that are two orders of magnitude smaller than those from the other two FNMs. Moreover, the mirror FNM solution exhibits a convergence rate that is increasing faster than the other two, as the horizon size decreases. The above two examples show that the mirror FNM works best at accurately enforcing local boundary conditions in PD models, especially for problems with singularities (in local models) along the boundary. In the following section, we will test the autonomous algorithm developed for the mirror FNM for problems with more complex geometries. 4.3. Steady-state diffusion in disks with single and multiple pre-cracks In this subsection, we will apply the autonomous mirror FNM to solve the PD formulation for diffusion in disks with cracks under Dirichlet boundary conditions. Using the algorithm developed in Section 3.4, the mirror FNM can be easily implemented for curved boundaries and cracks. 4.3.1. Disk with a single crack We first consider a disk with a single crack and boundary conditions imposed as shown in Fig. 19. The mirror node for each fictitious node, required by the mirror FNM, is determined by the algorithm described in Section 3.4. See Fig. 12 for how mirror nodes are found for those nodes with < 1 (which act as fictitious nodes) near the crack tip. Fig. 18. -convergence of PD results (using the different FNMs)to the classical value at point P (- 0.5, -0.48) in Fig. 15. 0 0.05 0.1 0.15 0.210-410-310-210-1100101102 Mirror FNM Taylor FNM Naive FNM 20 Fig. 19. Local boundary conditions imposed for diffusion in a disk with a pre-crack at the center. For verification, we also run a simulation of the corresponding PDE-based problem using the steady-state thermal solver in ANSYS workbench. Details of the ANSYS simulation, including how the “crack” is simulated, can be found in Appendix B. In the PD model, the crack is inserted by cutting all bonds that intersect with the pre-crack segment. The differences between the slightly different approaches as representing a crack between the two models should only have a trivial effect on the results. Contours of the results obtained by ANSYS and PD, with zoomed-ins around the crack region, are given in Fig. 20. The PD results match the FEM results closely, which shows that the autonomous mirror FNM works very well for problems with a curved boundary and with cracks. The jagged shape in the contours of the PD results appears because no smoothing interpolation technique is used in the visualization, we simply plot the obtained values at each PD node. 21 (a) (c) (b) (d) Fig. 20. Contours for the solution to problem shown in Fig. 19 obtained with: (a) ANSYSand (b) PD. In (c) and (d)we show zoomed-in regions around the crackfor the corresponding solutions. In practical simulations, the crack surface may have an angle with respect to the (uniform) discretization grid, whether it is a pre-crack or a new crack formed during the simulation. In order to demonstrate the generality of our algorithm for the mirror FNM, we solve the same problem but we rotate the uniform grid in the counterclockwise direction by 30° and 45° relative to the coordinates shown in Fig. 19, respectively, as shown in Fig. 21. Notice that when the crack is not aligned with the lattice, it may intersect with some PD nodes. These nodes are considered as free nodes (fully damaged) at which the constraint is the average of their family nodes (those not fully damaged). The PD solution with 30°-rotated and 45°-rotated grids are shown in Fig. 22. They match very well the solution obtained with the original lattice shown in Fig. 20. Therefore, the autonomous algorithm can deal with arbitrarily oriented uniform grids effectively. 22 (a) (b) Fig. 21. Part of the crack segment, shown by the red segment, and the PD grid after counterclockwise rotation by (a) 30° and (b) 45°. (a) (c) (b) (d) Fig. 22. Contours of PD solutions obtained with the mirror FNM for counterclockwise-rotated grids by (a) 30° and (b) 45°. In (c) and (d) we show zoomed-in views around the crack for the corresponding solutions. 23 4.3.2. Diffusion in a disk with crossing cracks A diffusion problem in a disk with two intersecting cracks is shown in Fig. 23. We use the same autonomous algorithm introduced in Section 3.4 without any changes to treat this case. A contour-plot for the PD solution and the zoomed-in picture around the two cracks with imposed Dirichlet BCs are shown in Fig. 24. Fig. 23. Diffusion in a disk with two intersecting cracks and associated local boundary conditions. (a) (b) Fig. 24. Contour-plot for the PD solution obtained with the mirror FNM (with the newalgorithm) (a) over the diskand (b) over the zoomed-in area around the intersecting cracks (crack lines are drawn only approximately on top of the plot). The examples shown in this section demonstrated that the autonomous algorithm for the mirror FNM works very well to enforce local BCs for complex geometries, including crack surfaces. Moreover, it can be readily employed for problems with moving boundaries such as corrosion damage, crack propagation, etc. 5. Conclusions We introduced an autonomous algorithm to find mirror nodes in the mirror-based fictitious nodes method (FNM) for arbitrary geometries, including domains with cracks. Without requiring an explicit mathematical description of the boundary, the algorithm is based on computing a nonlocal gradient, at fictitious nodes, to determine the “generalized” normal direction to the boundary for arbitrary geometry. The new algorithm allows us to easily generate the necessary data for mirror-based FNM to correctly impose the desired local 24 boundary conditions in peridynamic diffusion models, as well as reduce/eliminate the peridynamic surface effect caused by incomplete nonlocal region near the boundary. We compared the mirror-based FNM with the naïve and the Taylor-based FNMs for problems with or without singularities (in the corresponding local models) along the boundary, and showed that the peridynamic solution with the mirror-based FNM agrees with the classical solution best, especially for the problem with a singularity in the corresponding classical formulation. The other two methods showed “pollution” of the solution far from the location of the singularity. We applied the new algorithm to diffusion problems in domains with a curved boundary and with cracks. In these cases, too, the peridynamic solution obtained with our mirror-based FNM matched well the FE solution obtained in ANSYS. The same algorithm, with few modifications, should also work for other types of peridynamic models to enforce local- type boundary conditions, including those with moving boundaries and growing cracks. The new algorithm, which enables a fictitious node method imposition of local boundary conditions in peridynamic models on arbitrary domains, will allow for more accurate peridynamic solutions near the boundaries. High accuracy near arbitrarily-shaped boundaries and material interfaces is crucial in, for example, problems that involve crack initiation and propagation, or evolution of corrosion fronts. In such problems, the new algorithm introduced here will have a great impact. Acknowledgments This work has been supported in part by the US National Science Foundation CMMI CDS&E Grant No. 1953346 (program manager Joanne Culbertson) and by a Nebraska System Science award. This work was completed utilizing the Holland Computing Center of the University of Nebraska, which receives support from the Nebraska Research Initiative. Appendix A. Numerical implementation of peridynamic diffusion model with the fictitious nodes method For spatial discretization, we discretize the peridynamic body uniformly [26] into cells with nodes in the center of those cells. Fig. A.1 shows a 2D uniform discretization with grid spacing around a node . Non-uniform grids are also possible [20,53,54], which may conform better for shapes with, for example, rounded boundaries [49], but this is not pursued in this work. To discretize the peridynamic integro-differential equations, we use a meshfree method with one-point Gaussian quadrature [26] for the approximation of the integral term. Faster numerical methods such as the boundary-adapted spectral method [17][55] can be alternative options. 25 Fig. A.1. Uniform discretization for a peridynamic model. The circular region is the horizon region of node . The discretized PD Laplace’s equation (Eq. (4)) becomes: = 0 (A.1) The superscript means load step. The subscripts and denote the current node and its family node respectively, in the discretized domain. is the horizon region of node , includes all the nodes covered by (fully or partially), = , is the area of node covered by . The discretized versions for other equations are similar to Eq. (A.1) and are not given here. In Taylor-based and mirror-based fictitious nodes methods, the equilibrium system can be solved iteratively using the linear Conjugate Gradient (CG) solver combined with additional criteria to check for the convergence of the solution. At each iteration or solution step, the CG solver is called and the solution in the domain and fictitious region is updated, respectively. For the purpose of minimizing the overall computational cost, the tolerance in the CG solver is set to be 1e-2 at first and then decreased with solution steps by a factor until it reaches 1e-6. This treatment could make the simulation 50% more efficient than fixing the tolerance in the CG solver to be 1e-6 during the whole simulation. The system converges when the solution in the domain obtained between two sequential solution steps differ, in terms of norm-2 relative difference, by less than a given tolerance (1e-6 in this work). The detailed workflow for a complete simulation is shown in Fig. A.2. 26 Fig. A.2. Workflow for the peridynamic simulation with Taylor/mirror FNM. Appendix B. FEM modeling of the steady-state thermal problem To obtain the classical FEM-based solution for the problem shown in Fig. 19, ANSYS Workbench Steady- State Thermal solver is used. In the FE model, the two crack surfaces are generated by two arcs with the same small curvature and the maximum space between them is 0.01, which equals the grid size in the corresponding PD model. For the mesh, as shown in Fig. B.1, the element order is selected to be program controlled and all elements are triangles with the maximum size equals 0.05. The total number of nodes and elements are 9312 and 4548, respectively. All other options in the solver are set to be defaults. 27 (a) (b) Fig. B.1. FEM mesh (a) over the whole disk; (b) near the crack tip.