1. Introduction

Material properties emerge from phenomena on scales ranging from angstroms to meters. A multiscale treatment can provide a basis for an understanding of material behavior on different scales. The traditional multiscale approach couples two models operating at different scales. Alternative modeling strategy is called peridynamics, a novel non-local continuum model approach suggested by Silling, see ^[1], that models materials on a per-particle basis. This approach continualize the discrete models such as molecular dynamic models, by replacing inhomogeneities present on smaller length scales by an enhanced continuum description on larger length scales. Peridynamic models converge to the classical elastic material model, it can be used as an alternative approach for multiscale modeling as a single multiscale model valid over wide range of length scales, see ^[2] and ^[3]. In that way peridynamic can handle length scales beyond those typical for molecular dynamic models, opening for new possibilities for application to systems at mesoscale and macroscale. Later developments of the original theory, called state based theory of peridynamics, see ^[4], showed that any constitutive model from the classical continuum theory can be adapted to this state based theory of peridynamics.

The force interactions resulting from the discretized peridynamic equations are very similar to traditional molecular dynamic formulations, and this suggests that peridynamics can be implemented within a molecular dynamics framework code, see ^[5]. It was shown in ^[6] showed that peridynamics can be considered as an upscaling of molecular dynamics. The molecular dynamics freeware LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator), see ^[7], will be used for the simulations since this code also supports the peridynamics framework, see ^[8]. LAMMPS is a molecular dynamics simulator which is shipped with a peridynamics package PDLAMMPS, developed by Sandia corporation, allowing it to take a peridynamic approach to materials analysis. In addition to the bond based peridynamic model implemented as prototype microelastic brittle model (PMB), see ^[4] and ^[5], there are models for the state-based peridynamic linearly elastic (LPS), elastic-plastic (EPS) and viscoelastic (VES), see ^[5,9] and ^[10] solids implemented. Details about the implementation of the elastic-plastic and viscoelastic peridynamic models in PDLAMMPS can be found ^[11,12,13].

The peridynamic package in PDLAMMPS allows for modelling material restricted to isotropic behavior of material, limiting the application of the peridynamic theory. This creates a need for extending the existing models to include such behavior. An anisotropy extension of the elasto-plastic peridynamic model is presented in this work. The objective has been to develop anisotropic peridynamic model that easily can be implemented into PDLAMMPS. This was achieved by modifying the spherical influence function in the peridynamic models. A non-spherical influence function was introduced and implemented, i.e., the forces are not modulated based on the internal distance of the acting particles, but based on the direction of their relative position. This way, anisotropic materials can be modeled by aligning the model with coordinate system and thereafter instruct the influence function to modulate the forces differently depending on what direction they have relative to this coordinate system.

The model is introduced into source files of LAMMPS itself, extending the program in a straight-forward manner. The influence function being implemented takes internal forces of the object and truncates their magnitudes according to their direction in the coordinate system. The capabilities to represent directional dependent behavior are shown through simulating beams loaded in compression. The model is shown to be internally consistent, i.e., returning the the same results for the elastic-plastic (EPS) model as a spacial case of the newly implemented anisotropic elastic-plastic model (AEP). Clear qualitative differences in behavior is shown independently of pre-existing simulation parameters.

---

2. Theory of peridynamics

In the peridynamic theory the interactions between material points are expressed in terms of bond forces, permitting interaction between particles at a distance, in contrast to the classical continuum mechanics, where interactions between material points are expressed in terms of traction vectors, i.e., we assume a local interaction. Peridynamics treats internal forces within a continuous solid as a network of pair interactions, similar to springs. Each material point interacts with its neighborhood within a sphere of radius $\delta$, called the horizon, that serves as an internal length scale in the model.

Let $\textbf{x}$ be the position vector of the studied point and $\textbf{x}'$ the position vector of a point in its neighborhood, see Figure 1. The relative position of the two particles in the reference configuration is denoted by $\xi$ and the relative displacement vector with $\eta$,

$\xi = \textbf{x}'-\textbf{x} \ \ \ , \ \ \ \eta = \textbf{u}'-\textbf{u}$

(1)

where $\textbf{u}'$ and $\textbf{u}$ are the displacement vectors of the two points. $H_{\chi}$ is a sphere whose center is the studied point $\textbf{x}$. The radius $\delta$ is the horizon of $\textbf{x}$ and its horizon is the distance limit across which a pair of material points can interact, and $\delta$ can be interpreted as the length scale in the model, defined in the reference configuration in contrast to the cut-off radius in molecular dynamics, which is defined in the deformed configuration. The point $\textbf{x}$ interacts directly with the points in the sphere. The interactions with other points are considered negligible. In state-based peridynamics the equation of motion of a material point at $\textbf{x}$ then reads

$\rho({\bf{x}})\ddot{\bf{u}}({\bf{x}}, t) = \int_{H_x} \left\lbrace {\bf{T}} \left[ {\bf{x}}, t \right] \langle {\bf{x}}'-{\bf{x}} \rangle - {\bf{T}} \left[{\bf{x}}', t \right] \langle {\bf{x}}-{\bf{x}}' \rangle \right\rbrace dV_{{\bf{x}}'} + {\bf{b}}({\bf{x}}, t)$

(2)

Here, ${\bf{T}}\left[{\bf{x}}, t \right]$ is a force vector state field, $\textbf{b}(\textbf{x}, t)$ is the external body force density field, $\rho(\textbf{x})$ is the mass density in the reference configuration and ${\bf{u}}({\bf{x}}, t)$ is the displacement field. To ensure consistency with basic physical principles, it will be required that ${\bf{T}}\left[{\bf{x}}, t \right]$ satisfy balance of linear momentum for any bounded body. A limiting process $\delta\rightarrow0$ produces a classical constitutive model for the Piola stress, see ^[14]. The force vector state field ${\bf{T}}\left[{\bf{x}}, t \right]$ contains entirety the constitutive model. The force vector state ${\bf{T}}\left[{\bf{x}}, t \right] = t {\bf{M}}$ can be interpreted in term of mapping the bond between two point to a force per volume, where $t$ is the force scalar state and ${\bf{M}}$ is the deformed state.

We introduce a nonnegative scalar state defined on $H_{\chi}$ by $x\langle\xi\rangle = \| \xi \|$. The elastic-plastic force scalar state $t$ is split into two parts, the volumetric and deviatoric scalar force states

$t = t_i + t_d = -\frac{3K \theta}{m}\omega x + \alpha\omega(e_d - e_{dp})$

(3)

where $K$ is the bulk modulus, $\alpha$ is material constants, $e_d, e_{pd}$ are the deviatoric and plastic extensions of the material, respective, whereas $\omega$ is the influence function. The weighted volume $m$ is defines as

$m({\bf{x}}) = \int_{H_{\chi}} \omega\langle\xi\rangle x\langle\xi\rangle x\langle\xi\rangle dV_{\xi}$

(4)

The dilation $\theta$ is

$\theta({\bf{x}}, t) = \frac{3}{m({\bf{x}})} \int_{H_{\chi}} \omega\langle\xi\rangle x\langle\xi\rangle e({\bf{x}}, t) dV_{\xi}$

(5)

$e({\bf{x}}, t)$ is extension scalar field defined as $e({\bf{x}}, t)\langle\xi\rangle = \|\xi+\eta\|-\|\xi\|$.

---

3. Anisotropic influence function

The influence function $\omega$ in the peridynamic theory is used to weight the contribution of all the bonds. It regulates how the particles interact with each other and was first presented as a boolean function returning a value of either $0$ or $1$ for each bond, modulating the forces by simply setting them to zero in certain instances and leaving them unchanged in others, and were first introduced in the modern state-based formulation of peridynamics by. Different material behavior can be modeled by letting the influence function modulate the forces acting on particles in the simulation based only on their distance from each other. Such an influence function is called spherical, as forces have the same magnitude in every direction from the particle. In ^[15], it was also demonstrated that influence functions can be used to modulate nonlocal effects within a peridynamic model independently of the peridynamic horizon. Details on spherical influence functions, can be found in ^[6] and ^[15].

In this work, a non-spherical influence function, see Figure 2, was implemented in a peridynamic model—the forces are not modulated based on the internal distance of the acting particles, but based on the direction of their relative position as seen relative to some coordinate system. This way, anisotropic materials can be modeled by aligning the model with coordinate system and thereafter let the influence function to modulate the forces differently depending on what direction they have relative to this coordinate system. This isotropic model is modified to include anisotropy by introducing the ellipsoidal influence function and the rest of the model is unchanged, i.e., the same flow rule as in the original models is used.

The influence function being implemented takes internal forces of the object and truncates their magnitudes according to their direction in the coordinate system as follows: if the force ${\bf{f}}$ is directed exactly along a predetermined preferred direction ${\bf{u}}$, then ${\bf{f}}$ will not be modulated at all. If the force is instead directed along a second direction ${\bf{v}}$, orthogonal to ${\bf{u}}$, the force vector will be multiplied by a predetermined scalar $a < 1$. If the force is directed along a third direction ${\bf{w}}$, orthogonal to both ${\bf{u}}, {\bf{v}}$, then the force vector will be multiplied by a second predetermined scalar $b < 1$. If ${\bf{f}}$ instead points in any other direction, then ${\bf{f}}$ will be modulated such that its magnitude will be the distance from the center of an ellipsoid to its edge, said ellipsoid being defined as being centered in the origin, having a major axis with length $\|{\bf{f}}\|$ pointing in the direction of the preferred direction ${\bf{u}}$ and semi-axes of length $a\cdot \|{\bf{f}}\|$, $b\cdot \|{\bf{f}}\|$. Finally, the semi-axes are rotated about the vector ${\bf{u}}$ by some angle $\theta$. Thus, the only calculation to consider is that of finding the distance from the origin of a cartesian 3-d coordinate system to any ellipsoid centered in said origin.

---

3.1. Algorithm Design

The approach taken is return the distance to the edge of an ellipsoid with the major semi-axis being the length of the incoming force vector and having a predefined direction, the other semi-axes having a lengths some fraction of the major semi-axis. The pre-defined direction of the major semi-axis we call the preferred direction, and should correspond to the preferred direction of deformation of the solid being modelled. The ellipsoids are defined by inputting the preferred direction along with real numbers $< 1$ for the two minor semi-axes and a rotation about its own major semi-axis defined in radians. The algorithm calculates extrinsic rotations to turn the ellipsoid so that its major semi-axis points in the preferred direction of the solid, and then rotates it intrinsically (that is, rotates it around its own major semi-axis) to match the input rotation. This results in a $3\times 3$ rotation matrix $M$. The LU factorization of this matrix is then calculated, and the results stored in a single $3\ \times\ 6$ matrix. It is this LU factorization that is used in the calculations once the simulation is run. The distance from the center of the ellipsis to its surface is then calculated as explained in the next section for each bond interaction in the simulation. As the ellipsoid is scaled such that the longest such distance is along its major semi-axis and is the exact length of the force vector ${\bf{f}}$ being subjected to the function, its is obvious that the distance $d \leq \|{\bf{f}}\|$ depends on the direction of $f$ relative to the pre-defined preferred direction of the material. This will lead to larger internal forces in the preferred direction of the material and thus smaller deformations and stresses, along with larger deformations and stresses in other directions. A material can have more than one preferred direction—the algorithm simulates this by storing several ellipsoids, running the scaling routine on each one and then returning the largest value. An even more useful algorithm would weight the ellipsoids somehow, and use a smarter selection of which ellipsoid should be used to modulate the force vector.

---

3.2. Mathematics

The mathematical problem to be solved can be stated as follows: given an arbitrarily rotated ellipsoid in euclidean 3d-space, with arbitrarily chosen semi-axes, find where an arbitrary vector $\textbf{v}$ intersects the surface of the ellipsoid if the length of $\textbf{v}$ is also the length of the major semi-axis of the ellipsoid, thus ensuring that it will always intersect the ellipsoid. The method to do this is comprises two steps. First change basis so that the semi-axes of the ellipsoid becomes parallel with the basis vectors of the new coordinate system. Then normalize ${\bf{v}}$ and then using the equation for an ellipsoid find a scalar $d$ such that $d\cdot {\bf{v}}$ lies on the surface of the ellipsoid. More detail regarding the mathematics can be found in ^[16].

The first step is accomplished as follows: let ${\bf{a}} = (x_a, y_a, z_a)$ be the vector that constitutes the major semi-axis of the ellipsoid. Now, calculate how to rotate the basis vector ${\bf{e}}_x$ so that $(c\cdot{\bf{e_x}}) = {\bf{a}}$ for an arbitrary scalar $c$:

$Y = \arcsin(\frac{y_a}{\|{\bf{a}}\|}) \ \ \ , \ \ \ P = \arcsin(\frac{z_a}{\|{\bf{a}}\|})$

(6)

where $P$ is the "pitch" and $Y$ is the "yaw" of ${\bf{a}}$, ie the angles of ${\bf{a}}$ with respect to the basis of the coordinate system ("pitch" is what angle an airplane flying along ${\bf{e}}_x$ would have to have with the ground and "yaw" what angle it would have to turn left-right to start flying along ${\bf{a}}$). The ellipsoid is also rotated along its own major semi-axis, that is, around ${\bf{a}}$—this angle we call $R$ for "roll". These are the so called Eulers angels. We can create a rotation matrix to perform a rotation corresponding to each of these intrinsic rotations. These we call $M_P, M_Y, M_R$. An arbitrary rotation $M$ may be written as a product of three simple rotations:

$M = M_Y \cdot M_P \cdot M_R$

(7)

where $M$ is then a rotation matrix such that $M \cdot {\bf{e}}_x = {\bf{a}}$. We have

$\begin{equation} M_Y = \left[ \begin{array}{ccc} \cos(Y) & -\sin(Y) & 0 \\ \sin(Y) & \cos(Y) & 0 \\ 0 & 0 & 1 \\ \end{array} \right] \end{equation}$

(8)

$\begin{equation} M_P = \left[ \begin{array}{ccc} \cos(P) & 0 & \sin(P) \\ 0 & 1 & 0 \\ -\sin(P) & 0 & \cos(P) \\ \end{array} \right] \end{equation}$

(9)

$\begin{equation} M_R = \left[ \begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos(R) & sin(R) \\ 0 & \sin(R) & \cos(R \\ \end{array} \right] \end{equation}$

(10)

To change ${\bf{v}}$ into the desired basis, we simply do ${\bf{v}}' = M^{-1}{\bf{v}}$. Observe that the rotation matrix only rotates vectors, it cannot scale them. Thus, $\det(M) = 1$.

For the second step of the calculation, we note the general equation for an origin-centred ellipsoid with its semi-axes on the basis vectors:

$\dfrac{x^2}{a^2} + \dfrac{y^2}{b^2} + \dfrac{z^2}{c^2} = 1$

(11)

where $(x, y, z)$ is some point on the surface of the ellipsoid with major semi-axis $a$ and the minor semi-axes $b$ and $c$. We have changed basis into one where our ellipsoid is oriented in this fashion and are looking for the point where ${\bf{v}}$ intersects the ellipsoid to calculate the distance $d$ to it from the origin. Let ${\bf{w}} = \frac{\bf{v}}{\|{\bf{v}}\|} = (w_x, w_y, w_z)$. ${\bf{w}}$ is a vector pointing in the direction of ${\bf{v}}$ with length 1, and the sought distance $d$ is a scalar such that $d \cdot {\bf{w}}$ is a point on the surface of the ellipsoid found by traveling the direction of ${\bf{v}}$, i.e., the intersection of ${\bf{v}}$ with the ellipsoid surface. Thus, $d$ must be some multiple of $\|{\bf{w}}\|$ that satisfies the ellipsoid equation, so that

$\dfrac{(d \cdot w_x)^2}{a^2} + \dfrac{(d \cdot w_y)^2}{b^2} + \dfrac{(d \cdot w_z)^2}{c^2} = 1$

(12)

From the last equation the dustance $d$ to the surface of the ellipsoid can be calculated as:

$d = \dfrac{1}{\sqrt{\dfrac{w_x^2}{a^2} + \dfrac{w_y^2}{b^2} + \dfrac{w_z^2}{c^2}}}$

(13)

---

3.3. Implementation procedure

A c-style struct is created from user input and saved, comprising:

1. The preferred direction ${\bf{a}} = (a_x, a_y, a_z)$ for the modeled material as a vector of doubles. This is the vector that points in the direction of the major semi-axis of the ellipsoid used in the calculations.

2. Two numbers $b$ and $c$ such that $\|{\bf{b}}\| = b\cdot \|{\bf{a}}\|$ and $\|{\bf{c}}\| = c\cdot \|{\bf{a}}\|$, if the other two semi-axes are thought to point along vectors ${\bf{b}}$ and ${\bf{c}}$, represented as doubles.

3. A number $r$ representing the intrinsic rotation of the ellipsoid in radians.

4. A double called "weight", intended to use if extending this modification to allow for more specific handling of several ellipsoids.

This struct is then used to calculate the Euler angles $Y, P, R$ as described in the section above. These angles are used to construct the rotation matrices $M_Y, M_P, M_R$. The total rotation matrix $M = M_Y\cdot M_P \cdot M_R$ is then constructed. The matrices $L, U$ such that $L\cdot U = M$ are calculated and saved as a vector of vectors of doubles in the struct. This struct is saved in a vector of similar structs. All of this is done upon initialization of the simulation.

For every time step in the simulation, the influence function is called with a force vector ${\bf{v}} = (x, v_y, v_z)$ as input. The vector of structs is iterated over, using each LU factorization to calculate the distance from the center of the ellipsoid it represents to its edge in the direction of ${\bf{v}}$. These distances are temporarily saved in a vector of doubles. The largest number in this vector is then found and returned. If no preferred directions are input the vector of structs will be empty, and the influence function will instead return $\sqrt{v_x^2 + v_y^2 + v_x^2} = \|{\bf{v}}\|$, indicating an isotropic material.

---

4. Results

To explore the performance of the anisotropic influence function being implemented is tested by conducting a compression test on a beam. A beam with a cross-section of $11\times11$ particles and a length of 100 particles was constructed in PDLAMMPS, aligned along the x axis. The length of the beam is 30 cm and the cross section is $3\times3$ cm. The beam was loaded in compression by prescribing the displacement at the two end surfaces. The ends of the beam were set to move toward the center of the beam at 0.2 cm/s. The tests were run for a variable amount of time steps, enough for the beam to buckle and crack. In all simulation the damage parameter was set to 0.1 and the horizon set to 1.0 cm.

As the beam is loaded during the simulation, the bonds between the particles may break when they are stretched beyond a given limit. The damage associated with each particle can be computed in LAMMPS and represented by taking value between 0 and 1. The upper limit correspond to a complete broken bond that can never recover. In the results presented below the damage parameter and the stress in axial direction are compared for different simulation.

In the initial test simulation the major axes and the semi-axes were set to be equal, i.e., the model should return the same result as the original isotropic influence function in PDLAMMPS. The obtained stresses were accurate to six digits of precision when comparing the two models, the anisotropic influence function with the initial isotropic influence function. The next step of testing the internal consistency is to achieve identical results when given different ellipsoids with the same proportions, i.e., doubling the major axis but proportioning the semi-axes in the same way. This is to ensure that the non-spherical influence function modulates the force depending only on it direction.

The initiation of the slip band of peridynamic nanobeam with the original isotropic model and anisotropic model with two preferable directions are seen in Figure 3. With the isotropic model necking band is in the vertical direction, whereas when introducing anisotropic model with preferable direction, this can be controlled in the desirable direction.

Figures 4–7 consists of a series of tests with ellipsoid influence functions, where the major axes of the ellipsoids point in different directions in the x–y-plane of the simulation. Figure 4 shows result for the damage parameter and the stress in the loading direction for simulations with spherical influence functions, to use as a baseline to inspect the ellipsoid influence functions against and result with an ellipsoid influence function, where the major axis of the ellipsoid is directed at a $45^{\circ}$ angle from the x-axis. With the isotropic model necking band is in the vertical direction, whereas when introducing anisotropic model with preferable direction, this can be controlled in the desirable direction. The beam acts qualitatively differently when there is an ellipsoid influence function present compared to the spherical case. With all the ellipsoid influence functions the beam bends and cracks in the x–y-plane rather than in the x–z-plane, most probably corresponding to the fact that all the ellipsoids are oriented in, and thus give rise to anisotropy in, said x–y-plane. The approach taken here is replace the sphere with an ellipsoid by setting the distance to the edge of the ellipsoid with the major semi-axis being the length of the incoming force vector and having a predefined direction, the other semi-axes having a lengths some fraction of the major semi-axis. The pre-defined direction of the major semi-axis correspond to the preferred direction of deformation of the solid being modeled. Since the preferable direction in Figure 4 is $45^{\circ}$ angle from the x-axis, the maximum stress in the x-direction is lower compared to the isotropic model.

Figure 5 shows results with an ellipsoid influence function, where the major axis of the ellipsoid orthogonal and parallel to the loading direction along x-axis. The occurring slip bands are orthogonal to the the loading direction leading to much higher damage value for the case when the major axis of the ellipsoid orthogonal to the loading direction along x-axis. Comparing the axial stress it is clearly seen that higher stress level is reached when the major axis of the ellipsoid is oriented along the loading direction. On Figure 6 results with the preferable direction of $45^{\circ}$ angle from the x-axis in the x–y plane is compared to simulation with ellipsoid directed at a $22.5^{\circ}$ angle from the x-axis in the x–z plane. The axial stress is as expected higher for the ellipsoid directed at a $22.5^{\circ}$ angle from the x-axis. The fact that the preferable direction lies in different effects the deformation pattern of the beam.

The implemented model is equipped with probabilities to include several preferable direction. The model is implemented in a such way that generally any number of ellipsoids oriented arbitrary can be included. Several influence functions oriented in a desirable direction allows this way to model material behavior with preferable direction. Finally a simulation with two concurrent ellipsoids pointing in different directions was conducted. This is illustrated in Figure 7. In the first simulations the model is equipped with two ellipsoidal influence function, one with the preferable direction of $45^{\circ}$ angle from the x-axis in the x–y and the second orthogonal to the x-axes. second simulations the model is equipped with one with the preferable direction of $45^{\circ}$ angle from the x-axis in the x–y and the second parallel to the x-axes. The orientation of the two ellipsoidal influence function can be seen on Figure 7. The deformation pattern is similar but it difference of the values of the axial stress is clearly seen on Figure 7. The algorithm simulates this by storing several ellipsoids, running the scaling routine on each one and then returning the largest value. This leads to higher axial stress for the case when the second ellipsoid is orthogonal to the loading direction. When putting two ellipsoids in the influence function, a larger ellipsoid dominates a smaller one. This is working as expected—the associated forces modulated by an ellipsoid that envelops another one will always be the largest, and as such will always be selected. However, having only a part of the smaller ellipsoid extend out of the larger one will give rise to measurable differences. With a weighted approach, rather than selecting the largest force always, this could probably be made to be much more convincing.

---

5. Conclusions

An anisotropic peridynamic model has been formulated and implemented in PDLAMMPS. The anisotropy have been introduced by a non-spherical influence function replacing the original spherical influence function. To explore the effects of anisotropic influence function within the peridynamic theory, simulation of a beam loaded on compression have been performed. The suggested non-spherical influence function is not restricted to a specific peridynamic model and can in general be applied to extend any isotropic peridynamic model.

The model shows good consistency when simulations with equal major axes and the semi-axes are performed, i.e., the model return the same result as the original isotropic influence function implemented in PDLAMMPS. The model also shows good internal consistency by reproducing results accurately when given different ellipsoids with the same proportions. Doubling the major axis but proportioning the semi-axes in the same way is accurate to six digits of precision, which is as high as the visualizing software will display.

The beam acts qualitatively and quantitatively differently when there is an ellipsoid influence function present compared to the spherical case. With all the ellipsoid influence functions the beam bends and cracks in the x–y-plane rather than in the x–z-plane, due to the fact that all the ellipsoids are oriented in, and thus give rise to anisotropy in, said x–y-plane. The different orientations of the ellipsoids also give rise to different orientation of the slip bands of the buckling beams. The axial stress is affected in the preferable direction of the ellipsoidal influence function.

Finally we have demonstrated that two or more preferable direction can be introduced in the model and the influence functions can be used to control that. In that way the model can easily be extended to include several preferable direction to model elastic-plastic deformation of anisotropic heterogeneous crystalline matter.

---

Conflict of interest

The authors declare no conflict of interest.

---