Homogenization of the 1D Peri-static/dynamic Bar with Triangular Micromodulus

In peridynamics, boundary effects generally appear due to nonlocality of interparticle forces; in particular, end effects are found in 1D bars. In a previous work by Eriksson and Stenström (J Peridyn Nonlocal Model 2(2):205–228, 2020), a simple method to remove end effects in certain types of 1D bars, or to homogenize such bars, was presented for bars with constant micromodulus. In this work, which is a continuation of Eriksson and Stenström (J Peridyn Nonlocal Model 2(2):205–228, 2020), the homogenizing procedure is applied to bars with a linear, or “triangular,” micromodulus. For the examples studied, common in practice, the linear elastic behavior of a homogenized bar, is identical to that of a corresponding classical continuum mechanics bar, independently of the interparticle force range and total number of material points of the bar.


Introduction
Peridynamic theory, introduced by Silling [2], is a nonlocal extension of solid mechanics, in which each material point is connected to its neighboring material points through pairwise forces acting inside a closed horizon. The formulation allows handling of long-range forces in general and unguided modeling of fractures in particular [2,3]. Peridynamics is based on integral equations that only involve the displacement field, thereby avoiding spatial derivatives that do not exist across discontinuities in classical continuum mechanics.
The pairwise forces of two material points are related to displacement through the stiffness of the bond between the points. The bond stiffness is computed under the assumption that a material point has a complete neighborhood. However, the neighborhood of a point near a boundary is incomplete, and using the complete neighborhood bond stiffness in a boundary region, results in a softer material and larger strains in such regions.
For a 1D peri-static/dynamic bar, this effect appears near the ends of the bar, but affects the entire bar. In general, a very large number of material points, in the order of 100-1000, is required to simulate a homogeneous bar, i.e., to achieve a displacement field that does not deviate notably from the analytical solution of a classical continuum mechanics bar.
In a previous work [1], we derived the parameters necessary to remove the end effects in 1D bars with a constant micromodulus c, or to homogenize such bars. The micromodulus is the elastic stiffness of the bond between two material points. In this work, the homogenizing procedure is applied to bars with a linear, or "triangular," micromodulus. We use the term homogenization in the sense of finding a peridynamic discretization that acts like a material with constant properties in the local theory, and not in terms of heterogeneous material.
The constant and the triangular micromodulus, which are common shapes of c, was derived by Bobaru et al. [4] for continuous bars, i.e., bars with an infinite number of points inside the horizon. The shape of the micromodulus is known to affect wave dispersion and convergence rates [4,5].
The 1D bar has been studied over the years as a natural branch of the peridynamic theory. Silling et al. [6] examined the dynamics of an infinitely long bar and identified features not present in the classical theory, such as wavelike oscillations spreading out from the loading region. The dynamics of peridynamic infinite bars of homogeneous and heterogeneous media have also been studied by, for example, Weckner and Abeyaratne [7] and Mikata [8].
Bobaru et al. [4] numerically investigated three kinds of displacement convergence, for a given type of bar loading, with different number of material points and horizon size. The authors also studied adaptive grid refinement away from the bar ends, to facilitate multiscale peridynamic modeling. They found that the relative error in relation to the classical continuum mechanics solution of deformation decreased as the number of material points in the bar increased.
Madenci and Oterkus [9] studied the longitudinal movement of the midpoint of a peridynamic bar, constrained at one end and momentarily loaded at the other, and discretized into 1000 points. The purpose of the study was to provide a benchmark example to be used in subsequent 2D and 3D modeling.
Seleson and Littlewood [10] utilized polynomial influence functions (similar to kernels) to modulate the strength of nonlocal interactions of a peristatic bar of 100 points, with displacement constraints at the bar ends. The numerical error of displacements was reduced by half by employing the influence functions. Chen et al. [11] achieved also significant improvements in the convergence of finite peridynamic bars with the corresponding classical solutions by modifying the kernel.
Aguiar et al. [12] investigated the behavior of the displacement field near the ends of a peridynamic bar and found that depending on the micromodulus, the displacement field could become unbounded near the bar ends.
Nishawala and Ostoja-Starzewski [13] suggested an inverse approach to the analytical linear elastic deformation of 1D bars, by first assuming a given deformation and then determining the loading required to obtain the given deformation, resulting in polynomial functions for element-wise applied body forces.
In view of 2D, Le and Bobaru [14] compared eight methods for surface correction suggested by various authors. The maximum displacement difference between the peridynamic solutions with surface correction and the classical solution varied from single to double digits in percentage. In general, these eight methods assumed a continuum and depended on geometrical calculations to compensate for the truncated nonlocal neighborhood.
Prakash [15] used a least squares approach to calibrate on a bond level in 3D. The study obtained exact continuum mechanics behavior for material points next to a surface. However, the improvements at edges and corners were smaller than for other existing methods, as seen in Le and Bobaru [14]. Even though the method can be applied in 1D, the direct calibration of the equilibrium equations presented in Eriksson and Stenström [1], and in this work, is simpler and more straightforward.
The paper is set up as follows. The next section briefly introduces peridynamic theory and considers the derivation of the effective micromodulus c for a bar with a finite number of points inside the horizon. We then derive the parameters necessary to adjust bond stiffness near the ends of the 1D bar. We show that the values and the number of stiffness correction parameters depend on the extent of a point's neighborhood. We also derive the adjustments of the related stiffness matrix to facilitate easy numerical implementation in peri-static/dynamic modeling, for both the stiffness matrix formulation and the most common implementation in practice, looping over equilibrium equations formulation. In the final section, the results are presented and discussed.

Bond-Based 1D Peridynamics
The peridynamic equation of motion of a material point at position x at time t is given as [2,3]: where is the domain of the body, u is the displacement vector field, ρ is the mass density and b is a prescribed body force field. f is the pairwise force function (a vector) per unit volume squared, denoting the force exerted by the material point at x on the material point at x. This interaction between pairs of material points is called bond. The integral is defined over a region H x , of radius δ, called the horizon (Fig. 1). H x can be seen as a sphere, disk or interval, for 3D, 2D and 1D models, respectively. The material points inside H x , except x itself, are called family members of x. After a horizon size has been chosen, the material body is discretized by choosing a relative grid density factor m = δ/Δ, where m is the number of point spacings within the horizon "radius," or size, and Δ is the constant material point spacing. Convergence studies may be performed to justify the selection of horizon and discretization. In 1D, the number of points inside the region H x is given by N H x = 1 + 2m. For 2D and 3D problems, the relative grid density factor m, should have a ratio of at least 3 [3,9] and, in many cases, 4 or higher [5,16,17] to provide grid independent crack growth patterns.
A material is called microelastic if the pairwise force f between material points is derivable from a micropotential ω [2]: where ξ = x −x is the relative position of two material points in the reference configuration, and η = u x , t − u (x, t) = u − u is the corresponding relative displacement in the deformed configuration. A linear microelastic material results in the micropotential where s is the relative elongation of a bond: Differentiation of (3) according to (2) yields: (5) where (ξ + η)/||ξ + η|| = e is a unit vector along a line through the two end points of a bond in the deformed configuration. As we assume that a material point x does not interact with material points outside its horizon, f = 0 for ||ξ || > δ. The kernel c(||ξ ||)/||ξ || in the integrand in (1), is common in peridynamic mechanical problems. Other kernels are possible, and the selection influences the nonlocality, convergence, and thus the discretization applied [11].
In 1D axially loaded straight bars, the vectors η and ξ are parallel, thus implying the force expression where η i denotes the component of the relative displacement in the direction of the undeformed bond and ξ /||ξ || = e. Note that (5a) is the linearized form of (5).
where c 1 is a constant, ξ = ||ξ || is the distance from a subject point and δ is the horizon. The micromodulus is conveniently obtained by equating the peridynamic strain energy density to the classical strain energy density at a point embedded within the bulk of a material body.
For 1D bodies, with the micromodulus assumed to be constant, c(ξ ) = c 1 = 2E/(Aδ 2 ), and with the micromodulus assumed to be triangular, c 1 = 6E/(Aδ 2 ), which have been obtained by Bobaru et al. [4] for a continuous bar, i.e., for an infinite number of points N H x inside the region H x . The constant micromodulus c is independent of N H x only in two special cases: a) for N H x infinite and b) for N H x finite, only in the case where the horizon divides an element in two equal parts and the element weight is halved [1]. This feature does however not hold for the triangular micromodulus, which is independent of N H x only if N H x is infinite, but is always a function of N H x when N H x is finite. The micromodulus for an increasing finite N H x tends toward the micromodulus for an infinite N H x . The micromodulus' dependence on N H x is derived below.
For the triangular micromodulus, the question is raised whether the horizon can be drawn somewhere between two points, so that the triangular micromodulus for an infinite N H x , c ∞ = c 1 = 6E/(Aδ 2 ), also is independent of a finite N H x , i.e., so that c ∞ applies independent of the horizon size.
The relative point values of the triangular micromodulus are i.e., the first point value, associated with the subject point, is unity and the last value, associated with the point at or nearest to the horizon, is β/(m + β) (Fig. 2). For 0 < β ≤ 1/2, the strain energy density is where iΔ corresponds to ξ and Δ to dξ . Expand the sum is obtained. (9) is calibrated by taking W m equal to the classical continuum mechanics strain energy density Substitution of (10) into (9), using (11), yields from which with (7) we obtain and further For c 1 (m) to be independent of m and equal to c ∞ , we must have which can be simplified to For 1/2 < β ≤ 1, the strain energy density is The coefficient c 1 (m) in (18) is calibrated as above. Substitution of (19) into (18), using (11), yield from which with (7) we obtain and further For c 1 (m) to be independent of m and equal to c ∞ , we must have Equation (17) can be simplified to It is readily seen that m ≤ 0 for 0 ≤ β ≤ 1/2 in (17) and for 1/2 ≤ β ≤ 1 in (25). Thus, there is no positive m for which (17) and (25) are satisfied for 0 ≤ β ≤ 1. In other words, for the triangular micromodulus, the coefficient c 1 (m) is always a function of the number of points within the horizon.
Below are considered the two most common cases in practice: the horizon a) divides an element into two equal parts and b) falls between two elements.
Taking β = 0 in (14), that is, a horizon divides an element into two equal parts, gives The function on the right-hand side is set to f 1 (m): which decreases with m toward unity (Table 1). Thus, if c ∞ is applied to a bar with a finite number of points in the region H x , then m must be at least 10 to reduce the error in the micromodulus to less than one percent and at least 32 to keep it smaller than tenth of a percent. Putting where is a user's maximum acceptable error, it is readily deduced that the condition (28) is fulfilled whenever Taking β = 1/2 in (14) or (22), that is, a horizon is drawn between two elements, gives The function on the right-hand side decreases with increasing m toward unity (Table 2). Thus, if c ∞ is applied to a bar with a finite number of points in the region H x , then m must be at least 5 to reduce the error in the micromodulus to less than one percent and at least 15 to keep it smaller than tenth of percent. The condition (28) is fulfilled whenever Thus, the case when the horizon is drawn between two elements is more favorable (as opposed to dividing an element into equal parts). Note that the conditions considering the minimum number of points/element inside the region H x of a bar apply only to a bar with triangular micromodulus.

Homogenization
The equation of equilibrium of a subject point/element in a 1D linearized peristatic bar, obtained from the basic peridynamic equation of motion [2], is in standard notation [4] where u i is the only component of the displacement vector u of the material point situated at x i .  Certain characteristic properties of the 1D peristatic bar are discussed in Eriksson & Stenström [1]. Like the previous case, a 1D bar is loaded in uniaxial tension by equal and opposite forces applied at the end points of the bar. The displacement of the bar is accordingly symmetric and, in particular, the displacement at the mid-section of the bar vanishes.
To illustrate our procedure, we start with the simplest case possible, namely a bar with a minimum horizon size and a minimum number of points, so that the entire bar is just covered by the elements inside the horizons. We then generalize the procedure to include greater horizons, bars with an arbitrary, including infinite, number of points N . Both bars with full weight and half weight end elements are considered.

Bond Stiffness Adjustment
For a bar with δ = 2Δ or 5 2 , direct inspection of a bar end reveals that the bond sets are different only in the first and last point interspaces of the bar (but they are equal between themselves). Accordingly, we assume an adjusted bond stiffness between points 1 and 2 (but not of the bond between points 1 and 3, as this choice would also affect the stiffness in the interspace between points 2 and 3) and likewise between points 4 and 5.
For a bar with δ = 3Δ or 7 2 Δ, the bond sets in the first two and the last two bar interspaces are different. In particular, the bond set in the first interspace is different from that in the second. We assume therefore an adjusted bond stiffness between points 1 and 2, and different adjustment between points 1 and 3, but no more. The same holds for the last two interspaces.
For a bar with δ = 4Δ or 9 2 Δ, the bond sets in the three first and last bar interspaces are different. We assume therefore an adjusted bond stiffness between points 1 and 2, and further different adjustments between points 1 and 3 and between points 1 and 4, but no more. The same holds for the last three interspaces.

Introductory Example
For a bar with horizon δ = 2Δ, a horizon covers three points i = 1, 2, 3, the subject point and two family points. The relative point values of the micromodulus are shown in Table 3. Table 3 Relative point values of The range of a bond force on either side of a subject point is just one point interspace. The bond forces are just local, i.e. they only extend to the nearest neighbor on either side of a subject point, but, as all element weights are not equal, the bar is not inherently homogeneous.
The first three equilibrium equations for a bar with the smallest possible number of points for the horizons to cover the entire bar, N H x = 5, are from (33) where A = V ip /Δ. The first factor 1/2 in the first term in (2 ) is the relative weight of the left end element of the bar. Other fractions in (1 ) to (3 ) are the relative values of the micromodulus, taken from Table 3. The equilibrium equations of points 4 and 5 are "mirror images" of (2 ) and (1 ) and omitted here.
To homogenize the bar, an adjusted bond stiffness is assumed between points 1 and 2 and between points 4 and 5. Because of load and displacement symmetry, it suffices to consider just one half of the bar. Noting in particular that bond stiffness is adjusted in just the first two equilibrium equations, we get from (1 ) and (2 ) where α ij are dimensionless stiffness adjustment parameters to be determined. In bondbased peridynamics α 21 = α 12 = α.
Due to symmetry, the displacement u 3 = 0 at the midpoint of the bar. To homogenize the deformation of the bar, we further require that u 1 = 2u 2 and u 4 = −u 2 (39a,b) i.e., the displacement of point 1 is twice that of point 2 with respect to the midpoint of the bar. Inserting (39a) in (38) and solving for α, we get The equilibrium equations of a bar are conveniently written on the matrix-vector form Su = −b. In standard arrangement, that is, with row number equal to equation number and column number equal to subject point number, a diagonal element in S is the displacement coefficient of the subject point in the corresponding equilibrium equation. The coefficients of the family members are off-diagonal elements on the same row, taken in order.
With Su = −4b/(c 1 A) in the present case, the 3 × 3 submatrix S u of S, from which S for both a non-homogenized as well as a homogenized bar can be obtained by simple extension, invoking (1 ), (2 ) and (3 ). Note that for a bar with half weight end elements S u is non-symmetric. S u for a homogenized bar is obtained by substituting α according to (40) and for a non-homogenized bar, or a bar "as such," by taking the homogenizing parameter to unity, i.e., in this case α = 1 in (41). Similar procedures apply to all S u given below. Factors b/(c 1 A) are derived in detail in Appendix 1.

Second Example
Next, a bar with δ = 3Δ and N H x = 7 is considered. A horizon covers four points i = 1, . . . , 4, the subject point and three family points. The relative point values of the micromodulus are shown in Table 4.
The equilibrium equations of interest are Here, the first factor 1/2 in the first term in (2 ) and (3 ) is the relative weight of the left end element of the bar. The other ratios are relative point values of c taken from Table 4 and displacement denominators are the distance in number of point spacings from the subject point.
Assume a bond stiffness αc between points 1 and 2 and βc between points 1 and 3 (and likewise at the other bar end). From (1 ), (2 ) and (3 ), we get The parameters α and β do not appear in any of the remaining equations of concern. At the midpoint of the bar u 4 = 0. Homogenization requires that u 1 = 3u 3 , u 2 = 2u 3 , u 5 = −u 3 and u 6 = −u 2 (49a-d)  (47) and (48) and solving for α and β, we get α = 3 and β = 2 (50a,b) With Su = −12b/(c 1 A) in this case, the non-symmetric 4 × 4 submatrix S u becomes ⎡ from which the submatrix S u for a homogenized as well as a non-homogenized bar can be readily obtained.

Third Example
Finally, a bar with δ = 4Δ and N H x = 9 is considered. A horizon covers five points i = 1, . . . , 5, the subject point and four family points. The relative point values of the micromodulus are shown in Table 5.
The equilibrium equations of interest in this case are Assume a bond stiffness αc between points 1 and 2, βc between points 1 and 3, and γ c between points 1 and 4 (and likewise at the other bar end). (1 )-(4 ) give As before, α, β and γ do not appear in any of the remaining equilibrium equations involved.

Generalization
The homogenizing parameters for a bar with half weight end elements and varying horizon size are collected in Table 6. From the calculated parameters for m = 2 to 4, the parameters for m = 5 and m are inferred.

First Example
For a bar with horizon δ = 2Δ, a horizon covers three points i = 0, 1, 2, the subject point and two family points. The relative point values of the micromodulus (Table 3) show that the range of a bond on either side of a subject point equals just one point interspace. As the bond forces only extend to the nearest neighbor, of equal weight, on either side of a subject point, the bar is not nonlocal and therefore inherently homogeneous.

Second Example
Next, for a bar with δ = 3Δ and N H x = 7, the relative point values of the micromodulus are shown in Table 4. For m ≥ 3 the bond ranges are greater than one point interspace and such bars are thus nonlocal. Assume a bond stiffness αc between points 1 and 2 and βc between points 1 and 3 (and likewise at the other bar end). The bond stiffness is adjusted in just the first three equilibrium equations. From (42)-(45), in which the factor 1/2, accounting for the weight of a half weight end element, has been replaced by a factor unity to account for the weight of a full weight end element, we get the equations of interest (corresponding to (47) and (48)) The homogenization conditions (49a-c) result in α = 3/2 and β = 1 (67a,b) Note that for full weight end elements the stiffness matrix S is symmetric.

Generalization
The homogenizing parameters for a bar with full weight end elements and varying horizon size are collected in Table 7. From the calculated parameters for m = 3 and 4, the parameters for m = 5 and m are inferred. From the above examples, we conclude that for a bar with horizon m (point spacings), the number of solving parameters in homogenization is m − 1 and that the solving parameters for any m form the set Thus, the homogenizing parameters for a bar with full weight end elements is half the corresponding parameter for a bar with half weight end elements. This also holds for m = 2, as a unit parameter means no adjustment.

Matrix Construction
Inspection of the above equilibrium equations on matrix form, reveals that in the general case, the stiffness matrix S can be obtained from the unique (m + 1) × (m + 1) submatrix S u of S by simple diagonal expansion. For a bar with half weight end elements S u and S are non-symmetric, while symmetric for a bar with full weight end elements.
A general procedure to obtain the elements of the N × N stiffness matrix S for a bar of any number of points/elements, with half weight and full weight end elements, respectively, is described in Appendix 2.

First Example
For a bar with δ = 5 2 Δ and N H x = 5, the relative point values of the micromodulus are shown in Table 8.
The equilibrium equations are   Table 8 and displacement denominators are distance in number of point spacings from the subject point.

Second Example
Next, for a bar with δ = 7 2 Δ and N H x = 7, the relative point values of the micromodulus are shown in Table 9.
The equilibrium equations of interest are Assume a bond stiffness αc between points 1 and 2 and βc between points 1 and 3 (and likewise at the other bar end). The bond stiffness is adjusted in (1 )- (3 ). From (2 )-(3 ), we

Third Example
Finally, for a bar with δ = 9 2 Δ and N H x = 9, the relative point values of the micromodulus are shown in Table 10.
The equilibrium equations of interest in this case are Assume a bond stiffness αc between points 1 and 2, βc between points 1 and 3 and γ c between points 1 and 4 (and likewise at the other bar end).

Generalization
The homogenizing parameters for a bar with full weight end elements and different horizon size are collected in Table 11. From the calculated parameters for δ/Δ = 5/2 to 9/2, the parameters for δ/Δ = 11/2 are inferred.
From above examples, we conclude that for a bar with horizon m (point spacings) the number of solving parameters in homogenization is m − 1 and that the solving parameters for any m form the set

Homogenization of Bar with an Arbitrary Number of Points
An important result is that the sets of homogenization parameters (α), (α, β), (α, β, γ ), etc., according to horizon size, hold in general for a bar with an arbitrary number of points. To show this, consider a bar with an odd number of points and let n ≥ 0 be the number of points inserted between the midpoint of the bar and the m outermost points at a bar end (on either side of the midpoint), that is, N = 1 + 2(n + m). Then, for the midpoint Homogenizing requires proportional displacements satisfying all equilibrium equations.
To give an example, simple but general enough, insertion of (101) into (85) and (86) and canceling the common factor u n+3 yield 30α(n + 3) − (30α + 41)(n + 2) + 30(n + 1) + 9n + 2(n − 1) = 0 (102) 9β(n + 3) + 30(n + 2) − (9β + 71)(n + 1) + 30n + 9(n − 1) + 2(n − 2) = 0 (103) It is seen that factors of n sum up to zero in each equation and that the remaining terms yield the same solution as before α = 9/5 and β = 4/3 This essential result is verified by all cases considered above and therefore assumed general. It is left to the reader to verify that the general result holds also for a bar with an even number of points. The homogenizing parameters are thus independent of the number of points of a bar.

Results and Discussion
First, an overview of some characteristic results of this work is presented, followed by results and discussion of two case studies considered. The overview summarizes the deformation of the 1D bar types considered. The first case study concerns boundary/constraint conditions of a bar with full weight end elements, using the matrix solution, and the second case study considers a longitudinally vibrating bar with full weight end elements, using a quasi-static formulation by looping over the equilibrium equations. In both cases, the implications of different boundary conditions are considered.
The numerical implementation has been carried out in MATLAB. Static codes are inhouse and implemented as described in the previous section, and quasi-static codes are modified forms of readily available Fortran codes by Madenci and Oterkus [9].

Overview of Non-homogenized and Homogenized 1D Peridynamic Bars
In the overview, the total strain of entire non-homogenized and homogenized 1D bars, with half weight and full weight end elements and horizons sizes 2, 3 and 4 element interspaces, as well as 5/2, 7/2 and 9/2 interspaces, is compared with the corresponding classical continuum mechanics strain. The bars are loaded in uniaxial tension by equal and opposite forces, applied to the bar ends. The displacement of a bar, from which the total strain is obtained, is calculated with the static stiffness matrix solution method.
The total strain of a peridynamic bar in uniaxial tension is in general a function of the number of points in the bar, slowly decreasing with increasing number of points toward a limit for an infinite number of points, corresponding to the classical continuum mechanics strain. The strain in the bar also varies piecewise between points along the bar. These features are due to different bond sets between some few points at the bar ends, which number depends upon the horizon size. The end effects of a bar propagates through the entire bar regardless of the number of points in the bar [1]. In the homogenized bar, the end effects are removed and the deformation is equal to that of a classical continuum mechanics bar.  Figure 3 shows the results for three non-homogenized and three homogenized bars, with full weight end elements, δ = 4Δ, N H x = 9 and b/(c m=4 A) = (5/2)(σ Δ/E) for different boundary/constraint conditions. In Fig. 3a, the number of material points N a = 2N b − 1 = 17, in Fig. 3b N

Boundary/Constraint Conditions-Case 1
The homogenization parameters are α = 2, β = 3/2 and γ = 1. The Young's modulus of the bars is unity and they are loaded in uniaxial tension to unit stress. The length of the bars in Fig. 3b, c is half the length of the bars in Fig. 3a.
The bars in Fig. 3a are loaded at both ends by equal and opposite forces. This means that the displacement of the midpoint of a bar vanishes, i.e., u a 9 = u a (x 9 = 0) = 0, and that the deformation of the bar is symmetric with respect to the midpoint, i.e., u a (−x i ) = −u a (x i ).
The bars in Fig. 3b, c are loaded by a force at the right end. The bar in Fig. 3b is constrained at the left end by prescribing u b 1 = 0 and the bar in Fig. 3c by putting u c 1,2,3,4 = 0, i.e., by prescribing zero displacement over a length corresponding to one horizon radius at the bar end, here denoted Δ-constraint and δ-constraint, respectively.
For the non-homogenized bars, the displacements of overlapping parts are different in all three cases, i.e., u c (x i ) = u b (x i ) = u a (x i > 0), which means in particular that the bars in Fig. 3b and c are not symmetry halves of the bar in Fig. 3a, as in classical continuum mechanics. Thus, for a non-homogenized bar, traditional symmetry conditions do not apply. Displacement constraints applied in a boundary region inside the material body are suggested by Silling [2] and Macek and Silling [18] or in a fictitious boundary region of material points outside of the material body are suggested by Madenci and Oterkus and Bobaru et al. [9,19], with a typical length of Δ or δ [20]. Comparison of Fig. 3b and c shows that δ-constraint (Fig. 3c) results in displacements closer to the analytical solution than Δ-constraint (Fig. 3b).  For a homogenized bar on the other hand, u b (x i ) = u a (x i > 0), that is, the bar in Fig. 3b is indeed a symmetry half of the entire bar in Fig. 3a, but not the bar in Fig. 3c, for which u c (x i ) = u a (x i > 0). Thus, Δ-constraint allows symmetry properties of a homogenized bar to be exploited (but not δ-constraint).

Vibrating Bar-Case 2
A vibrating bar with a Young's modulus of 200 GPa is loaded by an initial strain condition of ∂u 0 /∂x = 0.001H (Δt − t), where H is the unit step function. The analytical solution of a bar constrained in one end is given by Rao [21]: where and ρ are strain and density, respectively. The peridynamic model parameters for the vibrating bar are the same as those for the static bar, i.e., δ = 4Δ and full weight end elements. Figure 4 shows the displacement of the element on the right-hand side of the constrained region. The constraint conditions of Fig. 4a-c correspond to the constraints in Fig. 3a-c. The bar length, discretization, constraint and the plotted element displacement are given in Table 13. It is seen that the stiffness of the vibrating bar and of the static bar are similar in Figs. 3 and 4a-c.
Lastly, we study the free end of the bar by plotting u 201 , u 100 and u 103 in Fig. 5. It is seen that the frequency of the vibration deviates for the non-homogenized solution in Fig. 5b. This deviation increases with increasing points within the horizon. To demonstrate this more clearly, δ is set to 6Δ in Fig. 5 instead of 4Δ. The time step size is set to 0.8 √ 2ρΔ/(2δAc m=4 ) (see Madenci and Oterkus [9] for details).

Conclusions
In this paper, the displacement and total strain of non-homogenized and homogenized 1D bars with triangular micromodulus and with half weight and full weight end elements, loaded in uniaxial tension, were derived analytically and computed numerically. The micromodulus is in general a function of the number of points inside the horizon both for bars with a constant micromodulus [1] and with a triangular micromodulus. For a given horizon, the dependence of the micromodulus on the number of points within the horizon is stronger for a triangular micromodulus than for a constant micromodulus. In the latter case, but not the former, the influence is very small and negligible in most cases in practice [1].
The total strain of a non-homogenized bar is in general a function of the number of points in the bar and of the number of points inside the horizon. The total strain slowly decreases toward a limit value with increasing number of points in the bar. Independent of the horizon size, this limit value equals the strain of a classical continuum mechanics bar, both for half weight and full weight end elements (f (∞) = 1; Table 12), contrasting the behavior of a bar with constant micromodulus for which the classical limit is obtained only for full weight end elements [1].
There are two steps in homogenizing a 1D bar with triangular micromodulus: (a) adjust the stiffness of the bar end bonds and (b) adjust the micromodulus, according to the number of points inside the horizon.
For a homogenized bar with triangular micromodulus, loaded in uniaxial tension, the total strain is always equal to the classical continuum mechanics strain, independent of the number of points in the bar and the number of points inside the horizon.
Classical continuum mechanics symmetry conditions do not apply to a non-homogenized bar. For half part of a symmetric model however, a symmetry boundary region of m elements/points with prescribed zero displacement is preferred over such a region of one single element/point (Figs. 3b, c, 4b, c and 5b, c), since the deviation from the full symmetric model is smaller.
On the other hand, half part of a homogenized symmetric model with a symmetry boundary region of one single element with prescribed zero displacement (Figs. 3b and 4b) gives exactly the same result as that obtained for a fully symmetric model, but does not so for a boundary region of m element/points. Funding Open access funding provided by Luleå University of Technology.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommonshorg/licenses/by/4.0/.

Appendix 1: Right-hand side b/(c 1 A) of bar end element equilibrium equations
For a horizon dividing an element into two equal parts we have from m = δ/Δ and (26) For a bar with half weight end elements b = 2σ/Δ on an end element. σ is applied stress, which yields and for a full weight end element b = σ/Δ, which yields For a horizon between two elements we have from (30) 1 and for a full weight end element The right-hand side of the bar end element equilibrium equations is written on the general form b 1 where C is a coefficient given in Table 14.

Appendix 2: Stiffness matrix
The equilibrium equations of the points/elements of a 1D bar with N elements are conveniently written on matrix-vector form Su = b, where S is a N × N stiffness matrix. This formulation allows the solution of a problem to be written formally as u = S −1 b, or in MATLAB code u = S\b. Due to conventional equation and point numbering, S presents a characteristic diagonal structure.
For a bar with full weight end elements, the first m rows are unique and the last m rows are "double reflections" of the m first rows. The rows in between are identical and show a repetitive pattern diagonally. The set of non-zero elements of row m + 1 is simply displaced one column to the right for every row-wise step downwards to row N − m. Downwards, a similar repetitive pattern is mirrored in rows N − m to N . The coefficients of row N − m + 1 are just copies of the coefficients of row N − m moved this time one step downwards and one step forwards (i.e., to the right) in matrix S. For rows N − m + 2 to N , after downwards-forwards diagonal copying, an outgoing coefficient, now at the right-hand side of the matrix S, is added to the diagonal coefficient of the row.

A2.1.2 Bar with half weight end elements
For a bar with half weight end elements, the procedure to determine the coefficients of the central equilibrium equations, with row number j = m + 1, . . . , N − m, is the same as above for full weight end elements.
The coefficients of row m are, in a first step, as above, upwards-backwards diagonal copies of the coefficients of row m+1. The coefficients in row m must however, in a second step, be adjusted due to the different weight of the left end element of the bar. Because an end element's relative weight is one half, the first coefficient in row m is halved and this value is added to the diagonal coefficient of the row.
Going upwards from row m, after upwards-backwards diagonal copying, the diagonal coefficient on row m − 1 is adjusted due to an outgoing element (as described above for full weight end elements). Further, due to the end element's half weight, the first coefficient in row m − 1 is halved and this value added to the diagonal coefficient of row m − 1. This procedure is repeated up to the second row. Thus, from row m − 1 and upwards, except for the first row, adjust the end coefficient Going downwards, a similar repetitive pattern is mirrored in rows N −m to N . Coefficient copying goes one step downwards and one step forwards in S. The last coefficient in a row is halved and the value is added to the diagonal coefficient of the row, except for the last row. From row N − m + 1 and onwards to row N , there is an outgoing coefficient at the right end of a row, and the diagonal coefficient is adjusted accordingly. The last coefficient in a row is halved and the diagonal coefficient is again adjusted.
Note, however, that the unique part of the stiffness matrix for a bar with half weight end points is not symmetric.
In summary, both for bars with full weight and half weight end elements, the stiffness matrix S of a bar, of any horizon size m and number of points N , can be constructed by expanding the following set of m+1 coefficients into a full row of coefficients, by reflecting