Two-dimensional implementation of the coarsening method for linear

: Peridynamic theory was introduced to overcome the limitations of classical continuum mechanics (CCM) in handling discontinuous material response. However, for certain problems, it is computationally expensive with respect to CCM based approaches. To reduce the computational time, a coarsening method was developed and its capabilities were demonstrated for one-dimensional structures by substituting a detailed model with a surrogate model with fewer degrees of freedom. The objective of this study is to extend the application of coarsening method for linear peridynamics for two-dimensional analysis. Moreover, the existing one-dimensional coarsening method was further explored by considering various different micromodulus functions. The numerical results demonstrated that coarsening approach has a potential to reduce the computational time with high accuracy for both one-dimensional and two-dimensional problems.


Introduction
Peridynamics is a nonlocal continuum theory developed by Silling [1] with the objective of overcoming the shortcomings associated with the Classical Continuum theory particularly as it relates to the ability to model fracture and damage in materials. To achieve this objective, the equation of motion of the classical theory was cast in an integral form. This paved the way for the unification of the mathematical modelling of continuous media, fracture and particles in a single modelling framework. Thus, peridynamics provided a robust structure that can find application in situations involving evolution and propagation of discontinuities such as crack nucleation and growth using the same field equation as in the continuous case. Another feature of the Peridynamic formulation is that it is a nonlocal theory that incorporates the concept of long-range forces by allowing interaction of particles located at a finite distance from each other through a pairwise force field.
The peridynamic theory has extensively been investigated and is shown to have capabilities of modelling a wide range of engineering problems [2][3][4][5][6][7][8][9][10]. However, for certain problems, it suffers from the disadvantage of being computationally more expensive compared to the classical continuum theory which is attributable to its nonlocal formulation. In an effort to mitigate against this issue, researchers have proposed a number of multiscale modelling techniques. An adaptive refinement method was proposed in [11] for bond-based peridynamic problems involving a variable horizon. A peridynamic-based framework for hierarchical multiscale modelling was proposed in [12] that couples molecular dynamics and peridynamics. These techniques allow coupling of models at different scales for effective inter-scale information transfer. A coarsening method was also proposed in [13]. The key idea of the coarsening method is to obtain accurate results by solving a reduced order model. The reduction in the order of model is achieved by substituting a detailed model with a surrogate model with fewer degrees of freedom. The objective of this study is to extend the application of a coarsening method for linear peridynamics proposed in [13] for two-dimensional analysis. Moreover, the existing one-dimensional coarsening method was further explored by considering different micromodulus functions other than that considered in [14].

Peridynamics
Recalling from the classical continuum theory, the equation of motion of a medium arising from conservation of momentum is as follows: (1) where is the mass density of the medium, is the stress tensor, is a vector of body force density, is the acceleration vector field of the material point within the medium at time . The key challenge in using Eq 1 to model discontinuous system behaviour is the presence of the divergence operator, which implies the existence of the spatial derivative of the stress field and consequently the displacement field within the domain of interest. However, since these field variables are not continuous over features such as crack tip and crack surfaces, the derivatives in such instance are undefined. In the peridynamic formulation, the equation of motion was casted such that the integral operators replaced the derivatives on the right hand side of Eq 1. A "bond-based" peridynamic equation of motion was originally proposed by Silling [1] as (2) where is the displacement vector field, is the neighbourhood of the particle located at point . is the pairwise force function with a unit of force per unit volume squared that a particle located at point exerts on a particle located at point .

Micromodulus function for linear peridynamics
If we assume a linear material behaviour, then the pairwise force function [1] takes the form (3) where is a tensor valued function called the micromodulus function given by (4) is the relative position of the particles in the reference configuration, and is their relative displacement. The peridynamic equation of motion (Eq 2) therefore specialises to (5) The micromodulus function is required to satisfy certain properties [1]. Moreover, the range of admissible functions that can be used is extensive [14]. However, only a few of them will be discussed here.
The first sets of examples of the micromodulus functions to be given belong to a class of materials called microelastic materials. A peridynamic material is said to be microelastic if the pairwise force function, , is derived from a scalar micropotential, : The micropotential, is the energy of a single bond and represents a local strain energy density [2] which is given by (7) where the stretch, s can be defined as (8) From Eqs 7 and 8, Eq 6 becomes (9) where is a function that is characteristic of the material and is called the bond constant. Differentiating Eq 9 with respect to and assuming small displacements yields (10) for two-dimensional configurations. For one-dimensional configurations, Eq 10 specialises to (11) In order to determine the bond constants for microelastic materials, the strain energy density of a given point can be defined [2]  In the case where the dependency of is in the form of an inverted triangular profile [3], in which case, bonds closer to the principal point are softer than bonds further away from it, such that, for bonds within the horizon. In this case, Eq 13 becomes

Two-dimensional micromodulus
In the two-dimensional case, let be the radius of sphere describing the horizon. Then, the two-dimensional equivalence of the micromodulus Eqs 17, 20 and 23 are respectively given as Notice that Eq 26 unlike its one-dimensional counterpart is not a constant function.

Coarsening method
The goal of the coarsening method is to provide a multiscale framework based on peridynamic theory that will predict macroscale response of a medium based on the underlying evolving microstructure, without necessarily resolving all microstructural details. This is achieved by a process of successive elimination of points from the medium called coarsening. Each successive coarsening results in a medium with a reduced material points as well as reduced level of geometrical detail. However, the properties of the coarsened medium are determined such that the effects of the excluded material points are implicitly included in the coarsened level simulation. The coarsening formulation presented here is adapted from [14] for the completeness of the paper.

Peridynamic coarsening formulation
To coarsen a detailed model as proposed in [13], let be a linear elastic peridynamic body as shown in Figure 2. Let be the set of linear admissible displacement fields on , and let be the micromodulus tensor associated with the material of the body where is the set of all second order tensors. Assume is a positive number that represents the maximum interaction distance for all points in , such that if Let , and let represent the set of admissible displacement fields on . and are called level 0 and level 1 body, respectively. The objective of this coarsening process is to express the internal forces acting on in terms of its own displacements only while implicitly accounting for forces points in exert on . To achieve this objective, let be an arbitrary point in and let be a positive number.

Discretization of the coarsening method
To carry out the numerical implementation of the coarsening method described in the preceding section, is discretised into nodes, which for simplicity are taken to have equal volume . Let be the position of node in level . For any nodes and , let The coarsened micromodulus is obtained by discretizing Eq 42 as To evaluate Eq 46, the resolvent kernel needs to be determined. A convenient means of obtaining the kernel function was proposed in [13]. Successive coarsening to higher levels can be achieved by following the same procedure as above, so that for any , To compute the coarsened micromodulus at level for node , the resolvent kernel must be determined. A method was suggested in [13] for determining the resolvent kernel. This method proceeds by creating and merging two sub-systems of equations into a global system of equations. To create such a system, two vectors representing the displacement field in the level body and representing the displacement field in level body are defined. The relationship between these two displacement fields is formulated as follows: The displacement fields and for nodes in level body are constrained to satisfy Eq 29: (48) The second system of equations is obtained by ensuring that the nodes complementary to nodes in level body satisfy Eqs 30, 31 and 32. Concatenation of these two systems of equations results in a global system of equations that has the following form in one dimension: where is the number of nodes in , is the number of nodes in and the diagonal elements are given by (50) In shorthand notation, Eq 49 can be written as where is square matrix in one-dimension and in two-dimensions. Let be the inverse of , in one-dimension will be the leftmost columns of . In two-dimensions, is the leftmost columns of .

Coarsening the micromodulus function
In this section, the coarsened form of various micromodulus functions will be presented. One-dimensional cases will be followed by a two dimensional case.

Coarsening of one-dimensional micromodulus functions
The one-dimensional problem is a homogeneous bar with a length of . The elastic modulus of the material is assumed to be . Three scenarios will be investigated by considering the following micromodulus functions: i. ii. iii.
The level 0 horizon is specified as . Note that this same horizon will be used in the coarsened models. The bar is discretized into 400 nodes with spacing Coarsening of the detailed model to the level body is carried out by retaining every fourth node of the level body. Similarly, coarsening of the level body to the level body proceeds by retaining every second node in the level body. Figure 3 shows

Coarsening of two-dimensional micromodulus function
The two-dimensional problem is a square plate with a thickness of . Let the elastic modulus of the plate to be and its micromodulus function is given by (55) Coarsening is carried out as shown schematically in Figure 7. Every second node in every second row of the level model is retained in the coarsened level model. Similarly, every second node in the second row of level is retained in level body. Figure 8 shows the micromodulus function of the level model as well as the coarsened micromodulus of the level model. It is noticed that the coarsened micromodulus is characterised by sharp peaks which is a reflection of the fact that it is only defined at the coarsened region. Figure 8a shows the response of the material point in the -direction for a displacement of the point in the -direction. Figure 8b shows the response of the material point in the -direction if the point is displaced in the -direction. Figure 8c shows the response in the -direction of the point when the point displaces in the -direction. Finally, Figure 8d shows the response of the point in the -direction when the point is displaced in the -direction.

Numerical results
This section comprises five numerical experiments to illustrate the application of the coarsening method in solving solid mechanics problems. The systems are assumed to be in a state of static equilibrium and thus an implicit static solution scheme has been utilized. The first three examples focus a one-dimensional bar under tension loading cases for homogeneous and composite materials with and/or without defect. The last two examples demonstrate the capability of the coarsening approach for two dimensional plate problems for isotropic and composite materials.

One-dimensional homogeneous bar under tension loading
This example is a numerical experiment aimed at illustrating the capability of the coarsening method in capturing the correct bulk behaviour of a system even with reduced degrees of freedom as opposed to if the degrees of freedom were to be reduced without the coarsening process. Consider a bar of length of material having a micromodulus function of the form: (56) The maximum interaction distance , where is the length of a material point. The detailed model (level 0) is discretised into points, so that . The elastic modulus of the material is . Coarsening of the level model to level proceeds by retaining every fourth node in level 0. Level 2 model retains every second node in level model as shown in Figure 3. A body force density of is applied to the rightmost level 2 node, while the third leftmost level 2 node is constrained to have zero displacement. The displacement field associated with the level 0 and coarsened levels 1 and 2 are shown in Figure 9. It is observed that even with reduced degrees of freedom as the detailed model is coarsened into level 1 and subsequently level 2, the displacement fields obtained from all these models are very similar to each other. This would not have been the case if the reduction in degrees of freedom were not achieved through the coarsening process.
The analytical solution based on the classical theory to the axial bar problem is of the form (57) where is the distance from the support, is the displacement at point , is the elastic modulus, is the cross sectional area, and is the force applied at the free end of the bar. The result from the analytical model and those from the detail level 0 model as well as the results from coarsened levels 1 and 2 model are almost identical.
To illustrate the need for the coarsening process in achieving model order reduction, the number of nodes was reduced from 1000 to 250 and further to 125 without following the coarsening procedure. As shown in Figure 9, although the results from the coarsened models and the results from the models with reduced number of nodes are very close to each other, coarsened model results agree better with the analytical solution.

One-dimensional composite bar under tension loading
In this example, a one-dimensional composite bar with a length of and a cross-sectional area of is considered. The composite bar has a periodic microstructure consisting of alternating strips and of the hard and the soft materials, respectively. Each strip is in length. The horizon for both hard and soft section is . The micromodulus of the composite system has the form: In other words, bonds with both ends in hard strip have hard properties while bonds with either ends in soft strip have soft properties. The elastic moduli of the hard and soft materials are and , respectively. The level horizon is specified as . Coarsening from the level to the level body is carried out as shown in Figure 3 by retaining every fourth node of the level body. Similarly, the level body is coarsened to level body by retaining every second node in the level body. A body force density of is applied to the three rightmost level nodes, while the three leftmost level nodes are constrained to have zero displacement. The computed displacements of coarsened level are shown in Figure 10. As expected, computation from the level revealed the most detailed microstructural information. As the computational regions get coarsened, these microstructural details get smoothed out. However, both detailed and coarsened models have similar global stretch. This demonstrates the fact that the effective properties produced by the coarsening process accurately reflect the bulk properties of the composite. where is a degrading term given by The horizon is specified as in the level and the bar is discretised into 200 nodes. The defect is located at the center of the bar as shown in Figure 11. The level contains every third node in the level and similarly, the level contains every third node of the level as shown in Figure 12. Prescribed displacement boundary conditions are applied to the three leftmost and three rightmost level nodes. The values of the prescribed displacements are given by , where is the coordinate of the node.  The displacement fields for the level model as well as the coarsened level and level models near the defect are shown in Figure 13. The three levels give identical displacement profile except for close to the site of the defect. This reflects the fact that due to fewer nodes in the coarsened models, wider spacing exits between nodes.  The displacement field obtained from solving both a detailed (level ) and the coarsened model (level ) are similar as can be seen from Figure 15 and Figure 16. The displacement profile of a specific grid of points along the -axis and -axis are plotted and shown in Figure 17 and Figure 18 respectively. It is observed that both models result in very similar displacement fields.

Two-dimensional composite plate under tension loading
This example will consider a two-dimensional square plate with periodic microstructure consisting of alternative strips and of hard and soft materials, respectively as shown in Figure 19. Each strip is in length. The horizon for both hard and soft sections is . The micromodulus of the composite system is given by (62) Figure 19. Discretisation and coarsening of the two-dimensional composite plate.
The horizon is specified as . The plate is discretized into 200 points in both and -directions with every point having . Coarsening of the detailed model and the applied boundary condition is same as that of previous example.
The computed displacement profile along the -direction for the identical boundary condition for the detailed model (level ) and the coarsened model (level ) is shown in Figure 20. It is observed that because the detailed model included more material points, it has been able to resolve more microstructural details than the coarsened models. However, both models have the same global stretch thereby demonstrating that the effective properties obtained using this coarsening method is an accurate representation of the bulk properties of the composite.

Computational cost
Determination of computational cost savings arising from coarsening a numerical model will follow from the idea suggested in [1]. Consider a body that is discretized into nodes in the level model. Let be the number of arithmetic operations used by a linear solver to solve the fully coarsened level model, where and are constants ( for Gaussian elimination [1]), and is the total number of nodes in level model. Suppose each level has as many nodes as the previous level, where is a constant, such that (63) where is the dimension of the problem (In the example on two dimensional homogeneous plate, and ). The computational effort in determining material properties at level from level is proportional to since in doing so, the inverse of the matrix has to be computed.
Therefore, for some positive constant , the total arithmetic operations that will be used to coarsen a numerical model and solve the fully coarsened level model can be written as It can be concluded from Eq 64 that the computational effort of coarsening up to level in a linear solver for a boundary value problem is reduced by a factor of over the computational effort needed to solve the detailed level model. For a two-dimensional problem such as the example on two-dimensional homogeneous plate, this factor will be . The computational price paid to determine the coarsened properties of the model is less than which is independent of .

Conclusions
The coarsening method developed in [13] extended for two-dimensional application in this work has demonstrated its robustness in modelling a range of problems. Moreover, several one-dimensional problems were considered to investigate the effect of different micromodulus functions other than considered in [14]. The examples solved in this work have demonstrated that this method is able to relate the bulk behaviour of materials to the details of their microstructure in two-dimensional problems. This allows us to reduce the order of the problem without losing accuracy in predicting the bulk response of the system. This translates to a reduction in computational time and storage required. Another advantage offered by the coarsening method is that the model retains its attributes such as boundary conditions and applied force density after the coarsening process. This is not the case if the model is reduced without going through the coarsening process. In this case, the boundary conditions and applied force density can no longer be identically applied as in the original model.
It is instructive to mention a number of caveats related to this method. This coarsening method is currently valid only for linear and static problems. This limitation is a consequence of assumption made in obtaining Eq 30.