A computational homogenization framework for non-ordinary state-based peridynamics

Peridynamic theory has been shown to possess the capabilities of describing phenomena that theories based on partial differential equations are not capable of describing. These phenomena include nonlocal interactions and presence of singularities in system responses. To exploit the capabilities offered by peridynamics in the homogenization of heterogenous media, a nonlocal computational homogenization theory based on peridynamic correspondence model (non-ordinary state-based peridynamics) is proposed. To set the development of the theory on a rigorous mathematical framework and to ensure consistency with the nonlocal nature of the peridynamic theory, a nonlocal vector calculus was used in the analysis of the nonlocal homogenization theory. The proposed theory is a two-scale micro–macro-homogenization strategy in which the constitutive relation at the macroscale is derived from explicit solution of a nonlocal volume constraint problem at the microscale. To justify the coupling between the two scales, nonlocal analogues of the stress and strain average theorems as well as the Hill–Mandel macrohomogeneity condition were derived. Validation of the proposed theory is achieved via numerical solution of Representative Volume Elements (RVE) from composite materials and comparing the results with those obtained by means of established methodologies.


Introduction
Understanding how materials behave has contributed tremendously to the development of modern technologies. Historically, the Cauchy continuum theory has been and continues to dominate as the tool used by scientists and engineers in studying the behaviour of materials. Central to the Cauchy continuum theory is the continuum hypothesis which presupposes that at an arbitrarily small scale (microscale), materials consist of a continuous distribution of infinitesimal particles that are idealised as point mass and each having physical properties such as mass, displacement and velocity which are volume averages over a finite-size domain defined in the microscale. The term microscale in this context should be understood to mean the scale at which the material is decomposed into continuum particles which does not necessarily have to be in the order of micrometres.
The state of the material at a point is governed by balance laws and interaction between the continuum particles is achieved through exchange of mass, momentum, and energy with immediate neighbours in accordance with Noll's principle of local action. The composition, geometries, distribution, and properties of these particles along with the characteristics of heterogeneities such as cracks and voids constitute the material's microstructure. The detail of this microstructure typically transcends many orders of magnitude. When developing a model to characterise these materials, a resolution threshold is selected below which information about the microstructure is ignored. If the selected level of resolution is not able to encapsulate important microstructural details, the model must be enriched to be able to capture the desired level of details.
A straightforward solution is the use of micromodels in which the resolution is refined until the model can explicitly resolve all important microstructural details. This solution strategy has been utilised to model microstructurally heterogeneous materials such as composites 1 3 using numerical techniques such as the Finite-Element method [1][2][3][4], Finite-Difference Method [5][6][7][8], and meshless methods such as the Element-free Galerkin method [9][10][11][12] to name just a few. This method of enriching the model suffers from several drawbacks amongst which include the fact that the microstructural details that plays important role in the response of the material may exist over wide orders of magnitude and explicit resolution of the microstructure for some applications may require computational resources that is prohibitively expensive. Other more fundamental problems include the fact that the response of a material specimen in the framework of the classical continuum theory is independent of the size and shape of the specimen [13]. This is, however, not the case as several studies [14][15][16][17] have indicated that material behaviour changes as their characteristic size decreases, a phenomenon known as size or scale effect [18]. Other problems include the presence of discontinuity in material response for which the classical theory is not good at handling [19] and resolution of nonlocal effects such as the existence of long-range interaction [20].
A variety of enrichment methodologies have been developed over the years to overcome the challenges associated with the classical Cauchy continuum theory. A broad class of these methodologies is called the generalised continuum theories [21] or microcontinuum theories [22]. In developing these methodologies, the local action hypothesis is relaxed or eliminated, and interaction between continuum points located at finite distance apart is permitted thereby leading essentially to continuum models that are nonlocal. The nonlocality implied here should be understood in the sense used in [23] to encompass 'weakly' nonlocal and 'strongly' nonlocal models. In this sense, the weakly nonlocal models, though mathematically adhering to the principle of local action, account for nonlocality by enriching the constitutive model with the first or higher gradient of relevant state variables or thermodynamic forces [24]. A foremost contribution to the development of the generalised continuum models is the work of Voigt [25]. Later contributions would lead to the development of a range of nonlocal models. The weak generalised continuum theories can broadly be classified into those that retain the kinematics of the classical Cauchy theory but account for nonlocality by incorporating higher gradients of the classical kinematic variables in the definition of the strain energy and those with extended kinematics, with the strain energy in this case expressed as a function of the classical kinematic variables and their gradient in addition to the newly incorporated variables. Among theories in the first category, we can cite the couple stress theories [26][27][28] and the strain gradient theories [29][30][31][32]. Of those in the second category, we can cite the Cosserat theory [33], the micromorphic theories [34], micropolar theory [35][36][37] and the internal variable models [38][39][40][41].
In the strong nonlocal category of enrichment, belong the integral-type nonlocal models. The development of these class of nonlocal models can be traced to the works of Kunin [42], Kroner [43] and Eringen [44]. In contrast to the weak nonlocal models, the integral-type nonlocal models completely jettison with the notion of locality even in their mathematical formulation by substituting the gradients of relevant state variables with a weighted integral of these state variables. By doing so, these models posit that the state variables at a continuum point is influenced directly by the state variables of all other points in the body or at least some finite neighbourhood of the point under consideration. Subsequent development of the integral-type nonlocal continuum theories includes the introduction of Peridynamic continuum theory [45], the Nonlocal Operator Method (NOM) originally proposed in [46] and later extended in [47] as Higher Order Nonlocal Operator Method (HONOM) to eliminate the requirement of shape functions in obtaining partial derivatives and a Long-range cohesive interaction model proposed in [48][49][50].
Peridynamics is a reformulation of the classical elasticity theory with extended capabilities of modelling discontinuous system response and long-range force interaction. Since its introduction, theory continues to receive growing interest from researchers and have thus found application in solving a wide range of engineering problems [51][52][53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69][70]. This is because it offers a theoretical and mathematical framework that allow for the modelling of systems with singularities without the need for modification of the governing equations. Discontinuous state variables are handled with the same equations as the continuous state variables. Another advantage offered by peridynamics is the introduction of an intrinsic length-scale parameter into the governing equation. This parameter allows the theory to be useful in modelling problems across a wide range of scales [71] and has been shown to correlate with the notion of an 'effective distance' introduced in [20] as a parameter that influences direction of crack in dynamic fracture.
Although peridynamics has been proven to overcome a lot of challenges associated with the classical Cauchy continuum theory, there are still problems for which resolution of relevant physics of the problem to obtain acceptable model fidelity require inordinate computational resources that is either prohibitively expensive or currently not available. The need to find a balance between acceptable model fidelity and lowering computational cost has motivated the development of a range of multiscale enrichment protocols for the peridynamic theory. The objective of these protocols is to enhance the capability of peridynamics in capturing and utilising information across spatial scales at acceptable computational cost. This objective is achieved in one of two approaches. The first approach is by developing a framework that allow the coupling of peridynamic theory with other theories of mechanics with the goal of leveraging on the strengths of the coupled modelling theories. This will fall under the purview of concurrent multiscale modelling. The second is achieved through frameworks that allow coupling of models at different length scale in the same part of a computational domain and are categorised as Hierarchical multiscale modelling.
A variety of concurrent multiscale modelling framework for peridynamics have been proposed. These approaches can broadly be categorised into monomodel and multimodel approaches. In the monomodel approach, peridynamics is used to model the entire computational domain. Multiscale capability is achieved by refining the grid appropriately to resolve details at different length scales. Development in this respect includes the adaptive grid refinement method proposed in [72][73][74]. Refinement of the grid requires changing the grid density and the horizon which in turn induces a change in the micromodulus function of the peridynamic model. This change is done relying on a scaling algorithm. A variable horizon method was proposed in [75] to solve the problem of ghost forces that is a consequence of using variable horizon by introducing the concept of partial stress and a splice technique. To resolve the problem of ghost forces without recourse to the partial stress and splice technique of [75], a dual horizon method [76] was proposed with extended capability of covering both bond-based peridynamics and state-based peridynamics. The dual horizon method was extended as Voronoi-based peridynamics in [77] for non-uniform discretization based on Voronoi diagrams.
The multimodel concurrent approaches are hybridization frameworks that focus primarily in developing coupling methodologies that allow peridynamics to be used concurrently with other modelling frameworks such as the classical continuum or atomistic modelling theories. The goal is to leverage on the advantages offered by each of the coupled models. These methods can broadly be classified as kinematic-based, force-based or energy-based coupling methods. In the kinematic-based methods [78][79][80][81], coupling is achieved by implementing a contact algorithm over the interface region that requires the satisfaction of a series of kinematical constraints. The force-based approaches [82][83][84][85][86][87][88] rely on the balance of force in the interface region to achieve coupling of the models. In some cases, nonlocal weight functions are used to account for the contributions from the models in the balance of force. In the energy-based approach, coupling relies on the principle of energy conservation in the interface region. Two notable subclasses of the energy-based approach have been proposed such as the morphing methods [89][90][91][92][93][94] and the Arlequin coupling method [95].
In the hierarchical multiscale subclass, efforts have been expended by several researchers that led to the development of different hierarchical methods. A class of the hierarchical methods that we will designate as Model Reduction techniques in the framework of peridynamics was first proposed in [96] as coarsening method. The key idea of the Model Reduction techniques is to be able to capture the behaviour high-fidelity models using fewer degrees of freedom. The reduction in the order of model is achieved by substituting a high-fidelity model with a surrogate model with fewer degrees of freedom. The coarsening method was extended in [97] for two-dimensional applications. A related method of model reduction based on static condensation of a highfidelity peridynamic model was proposed in [98].
Another subclass of the hierarchical methodologies developed for peridynamics is the homogenization methods which encompasses several related methods. The earliest attempt at developing a homogenization framework for peridynamics appeared in [99,100] in which a two-scale solution expansion of the peridynamic equation for a heterogeneous medium was proposed based on the concept of two-scale convergence. The framework is a three-step solution strategy of computing an average displacement field as a solution of a peridynamic macroscopic equation and computing micro-level displacement field from a solution of a microscopic equation. The displacement field of the heterogeneous medium is found by superimposition of the micro-displacement field unto the macro-field in a final step. Although the mathematical framework was developed, numerical validation for the method is yet to be done.
A class of mean-field homogenization methods that can be categorised as a generalisation of the self-consistent methods [101] of the classical theory to peridynamics were proposed in [102][103][104] by extending the effective field hypothesis of the classical elasticity of composites to the nonlocal peridynamic framework. The effective elastic properties of a heterogeneous medium are determined through the introduction of a stress polarisation tensor.
Among the homogenization methodologies exists a class designated as full-field (computational) homogenization techniques which achieve a much higher resolution of the microscopic fields than the mean-field methods. In several works [105][106][107], various computational homogenization methodologies have been proposed to evaluate effective properties of heterogenous medium in the bond-based peridynamic framework based on the description of the microgeometry of the heterogeneous medium by solving a well-defined microscale volume-constrained problem. An extension of the computational homogenization framework to Ordinary state-based peridynamics was proposed in [108].
The objective in this study is to develop a multiscale constitutive theory for non-ordinary state-based peridynamics (peridynamic correspondence model) in the framework of computational homogenization. The methodology will be set on a rigorous mathematical foundation in a manner that will account for the nonlocal nature of the peridynamic theory.

3
To achieve these goals and for the sake of self-contained presentation, this paper will start off with a brief review of relevant results from a nonlocal vector calculus developed in [109]. This provides a mathematical framework that is consistent with the nonlocal nature of Peridynamics. Key results from the nonlocal vector calculus will be used as the building block to derive the nonlocal Peridynamic model in Sect. 3. Discussion on the discretization methodology for the peridynamic model adopted in this presentation is provided in Sect. 3.3. The development of a nonlocal computational homogenization for the peridynamic model is presented in Sect. 4 while numerical experiments to validate the proposed methodology is presented in Sect. 5. Concluding remarks and discussion on future research outlook is presented in Sect. 6.

Nonlocal vector calculus primer
In developing the nonlocal vector calculus, two types of functions and operators were defined in [109]. Let k, m and n be positive integers and let x and x ′ be points in ℝ n . For a given domain Ω ⊆ ℝ n , functions or operators that maps Ω into ℝ m×n or ℝ n or ℝ are called point functions or operators, respectively. On the other hand, functions, or operators from Ω × Ω into ℝ m×n or ℝ n or ℝ are called two-point functions or operators, respectively. Point functions and two-point functions could be scalar, vector or tensor valued functions.
A very important concept to start with in this review is the nonlocal flux. Given a tensor two-point function ∶ ℝ n × ℝ n → ℝ k , then the definition is the nonlocal flux of q from Ω 1 into Ω 2 where ∫ Ω 2 x, x � dx � is identified as the nonlocal flux density into the region Ω 2 from point x ∈ Ω 1 . It can be deduced from (1) that the nonlocal flux is not necessarily zero even if the intersection of the closures of Ω 1 and Ω 2 is an empty set. This is in stark contrast with the local flux which is zero if Ω 1 ∩ Ω 2 = ∅ . The nonlocal flux density is related to the intensive quantity q through a constitutive relation. If x, x ′ is assumed to be antisymmetric, then the following statements are true: 1. There is no self-interaction, i.e.
2. The nonlocal action-reaction principle holds for The action-reaction principle given by (3) simply states that the flux from Ω 1 into Ω 2 is equal to the flux that exits Ω 2 into Ω 1 .

Nonlocal divergence and gradient operators and their adjoint
Given the two-point function v ∶ ℝ n × ℝ n → ℝ k and the scalar two-point functionu ∶ ℝ n → ℝ . Let x, x � ∶ ℝ n × ℝ n → ℝ m be an antisymmetric vector twopoint function. The action of nonlocal divergence operator D and its adjoint D * on v and u , respectively, are defined as Given the scalar two-point function ∶ ℝ n × ℝ n → ℝ and the vector point function u ∶ ℝ n → ℝ k , for a given antisymmetric vector two-point function x, x � ∶ ℝ n × ℝ n → ℝ m , the action of the nonlocal gradient operator G and its adjoint G * on and u, respectively, are defined as Observe that, unlike in local calculus which deals with point functions only, nonlocal calculus involves two kinds of functions: point and two-point functions. This, therefore, necessitates the definition of alternative forms of the nonlocal operators defined in (7) and (8). The alternative forms of the nonlocal divergence and gradient operators were given in [28]  It is possible to apply the nonlocal divergence operator on a tensor function and the nonlocal gradient operator on a vector function. Let x, x � ∶ ℝ n × ℝ n → ℝ m be an antisymmetric vector two-point function. Given the tensor two-point function ∶ ℝ n × ℝ n → ℝ m×k and the vector twopoint function v ∶ ℝ n × ℝ n → ℝ k , the nonlocal divergence is defined by its action on as The action of the nonlocal adjoint operator D * on v is given by The action of the nonlocal adjoint operator G * on is given by

Interaction kernels and domains
In (7) and (8), the two-point vector functions x, x ′ and x, x ′ are also known as the interaction kernels. In the context of Peridynamics, these kernels are assumed to have a finite domain that does not map to zero. Given two points x, x � ∈ ℝ n , and ∈ ℝ + , let B (x) be a ball or radius centred at x , then for example, the interaction kernel (x, where is the interaction radius also known as the horizon in the context of PD. Interaction kernels that satisfy (12) are called truncated kernels [110] or localised kernels [111]. Another key concept that is connected to the notion of truncated kernels is the interaction domain. Let Ω ⊂ ℝ n be a bounded open set. The interaction domain Ω I consists of points outside of Ω that interact with points in Ω . Given the truncated interaction kernel (x, x � ) , an interaction domain Ω I is defined as The interaction subdomain contains all the points x ′ in the complement domain ℝ n �Ω that interact with points x in Ω . Many geometrical relationships exist between Ω and Ω I [109]. A typical such relationship is shown in Fig. 1.

Nonlocal interaction operators
Given a domain Ω , let Ω I be the interaction domain associated with Ω as defined in Sect. 2.2. Corresponding to the nonlo- Corresponding to the nonlocal gradient operator ( )(x) , a point interaction operator S( )(x) ∶ ℝ n × ℝ n → ℝ m×k is defined as

Nonlocal integral theorem
A very important outcome of the nonlocal operators developed in the proceeding subsections is statement of the nonlocal Gauss theorem. Recall from (2) that

Fig. 1 Interaction domain
Consider the right-hand side of (16). If the kernel function x, x ′ is antisymmetric and (2) holds, then: From (18) and considering (3), it can be deduced that Thus, (17) is the mathematical statement of the nonlocal Gauss theorem which postulates that the integral of the nonlocal divergence of v over Ω is equal to the total flux exiting Ω into Ω I . We next consider the nonlocal analogue of integration by parts expressions involving the nonlocal divergence and gradient operators. Given the point functions

Nonlocal weighted operators
Since the classical continuum theory is a local theory, the functions utilised are point functions. This is not the case with peridynamics which is a nonlocal theory, hence some of the functions utilised are point functions while some are two-point functions. The differential operators defined in Sect. 2.1 are operators that act on two-point functions. To complete the definition of the nonlocal operators, another class of operators need to be defined that act on point functions. Given the point function u ∶ ℝ n → ℝ k and the point functionv ∶ ℝ n → ℝ , let x, x � ∶ ℝ n × ℝ n → ℝ + and let the operators D and G be as defined in Sect. 2.1, then the weighted nonlocal divergence operator acting on u is defined as The weighted nonlocal gradient operator acting on v is defined as The extended application of the weighted nonlocal divergence and gradient operators on tensor and vector fields, respectively, follows as in (8) and (10). In addition, as is with the unweighted nonlocal operators, adjoint operators can also be defined for the weighted nonlocal operators. Let D * and G * be as defined in Sect. 2.1 and ∶ ℝ n × ℝ n → ℝ be a non-negative symmetric function known as the weight function, then the action of the adjoint operator D * corresponding to the nonlocal weighted divergence operator on v is defined as The action of the adjoint operator G * corresponding to the nonlocal weighted gradient operator G * on v is defined as The relationship established between D and −G * and between G and −D * allows us to, respectively, write (21) and (22) as and Equations (25) and (26) serve to define two forms to each of the weighted nonlocal divergence and gradient operators, respectively, and have been shown to be equal [109]. For example, from (25), the first form of the weighted nonlocal divergence is given by and the second form is given by

Nonlocal differential operators for peridynamic application
Often, the governing equations in the classical continuum theory are expression of physical balance laws that are composed of partial differential operators. Although the peridynamic theory as a nonlocal model replaces the spatial derivatives with integral operators, it can be shown to have retained the structure of the balance law in the classical theory. This is achieved by introducing nonlocal analogues of the differential operators used in the local theory. The nonlocal operators presented in Sect. 2 will be used to derive application specific to nonlocal differential operators and the generalised notion of nonlocal derivative for application in peridynamics.
then the nonlocal unweighted and weighted divergence and gradient operators from (4), (6), (25) and (26) are, respectively, given as where (32) and (33) appeared in [113] are defined as the notions of nonlocal material divergence and gradient operators, respectively. From (28), the second form of the weighted nonlocal divergence operator is

Peridynamic model
Peridynamics is a nonlocal alternative to the classical continuum mechanics developed in [45]. The objective was to create a nonlocal continuum theory that is capable of modelling both continuous and discontinuous material response in a single framework. This was achieved by replacing the differential operators in the equilibrium equation of motion of the classical continuum theory with an integral operator. Given a bounded open domain Ω ∈ ℝ n , in peridynamics, a continuum point x ∈ Ω interacts with infinitely many other points located within its domain of influence. If this domain of influence is assumed to be a ball B (x) of radius > 0 centred at x , then is called the horizon of x , such that where B (x) is called the family of x . Interaction between two points x and x ′ is called a bond and the distance | | = |x � − x| in the undeformed reference configuration is called the bond length. The relative displacement in the deformed configuration is given by where u x ′ , t and u(x, t) are the displacements of material points at x ′ and x, respectively, at time t.
In the original development [45] that came to be known as the Bond-Based Peridynamics, the force in each bond is a central force and each bond is independent of all other bonds. This lack of dependence limits the value of Poisson's ratio to 1/3 for 2D and 1/4 for 3D isotropic solids. To circumvent this limitation and obtain a more general material model, State Based Peridynamics was developed. To pursue the development of the state-based peridynamics, [112] Perhaps the most important state in the peridynamic formulation is the vector state and there are three important vector states that are worth mentioning here: the reference position vector state, deformation state where y ′ x ′ , t and y(x, t) are, respectively, the coordinates of the points x ′ and x in the deformed configuration at time t . The force state _ is a function that is associated with each bond in the family of point x with some force density vector such that A very import tensor (already introduced in Sect. 2.6) in the formulation of the state-based peridynamics is the shape tensor . This tensor can also be defined using the notion of states. Following from (37), the shape tensor is defined as where _ is as defined in (38). Although vector state and second-order tensors both map vectors to vectors, they are essentially different. For example, a state is in general nonlinear in its argument and is infinite dimensional. In contrast, a second-order tensor is linear function of its argument and has a dimension of 9. It was, however, demonstrated in [31] that given a second-order tensor , it is possible to obtain a vector state through an expansion operation defined as where ϵ _ ( )⟨ ⟩ is the vector state expanded from . Conversely, given a vector state _ ∈ V , a second-order tensor can be obtained by a reduction operation defined as The governing equation of motion in PD can be formulated using the nonlocal vector calculus presented in Sect. 2 based on a statement of balance law which postulates the dependence of the rate of change in the content of an extensive quantity over a given domain on the rate at which the quantity is produced within the domain and a flux through the boundary of the domain.  (13), then: Applying the nonlocal Gauss theorem (17), then (46) can be written as Thus, (44) becomes From localization theorem, (48) can be written as From (34), (49) can be written as where (50) follows from the antisymmetric property of . From (42), the integrand of the integration in (50) can be understood to be a vector obtained when a vector state ϵ _ K −1 is expanded from a second-order tensor K −1 acting upon the bond xx ′ or x ′ x as the case maybe. If we write (50) can be written as Equation (52) is the state-based peridynamic equation of motion.
_ [x, t] and _ x ′ , t are the force states at x and x ′ , respectively. When these states act upon the bonds xx ′ and x ′ x , respectively, the results are bond force density vectors acting at points x and x ′ , respectively. Notice that we still do not know what the second-order tensor represents. Its identity can be established from a constitutive material model that will relate the bond force density to the displacement. This will be the subject of Sect. 3.2.

Nonlocal kinematic quantities
In this section, the nonlocal differential operators in Sect. 2.6 will be used to derive the expression of relevant kinematic objects that are nonlocal analogue of their counterparts in the local theory.

Gradient of the displacement and deformation vector fields
An important quantity in the formulation of continuum mechanics is the gradient of the displacement vector. Let u(x, t) be the displacement of a point x at time t . Then, from (33), the gradient of the displacement field at x as a function of the undeformed bond is given as Another key kinematic quantity that plays an important role in the development of the state-based peridynamic theory is the concept of the deformation tensor. In the nonlocal peridynamic setting, the deformation tensor, denoted as, F, is the finite-dimensional equivalent of the deformation state. Let L + be the set of all second-order tensors with positive determinants. Let Ω 0 and Ω t be the reference and deformed configuration of a body undergoing deformation. Let x be the position of a material point in Ω 0 and y be its position inΩ t . Let ∈ L + exists such that the deformed image of the bond xx ′ is given by Using (33), the nonlocal material gradient of the deformation as a function of the undeformed bond is evaluated as Equation (55) is the definition of the nonlocal deformation gradient given in [112]. We can further write (55) as where the second equality follows from the linearity of the integral operator. Notice that the nonlocal deformation gradient is devoid of any notion of the local derivative operator which would have required that the deformation field be continuously differentiable (at least in a weak sense for the case of Finite-Element Method). In its present form, the deformation gradient is still defined in the presence of singularities such as cracks.

Nonlocal strain tensor
We can define the notion of strain by comparing the length of a bond before and after deformation. Let = x � − x and = y x � − y(x) , then from (38) and (54), we have If we define a deformation or strain matrix as t).
then (57) can be written as where is the nonlocal analogue of the Green-Lagrange strain tensor and is the nonlocal deformation gradient defined in (55). Considering (56), (58) can be written as Following from the assumption of infinitesimal deformation where the displacement gradient is small, that is | ≪ 1 , we may neglect the nonlinear term in the definition of the Green-Lagrange strain tensor so that (60) reduces to where is the infinitesimal strain tensor and G s denotes a symmetric tensor operator.

Constitutive model
We will now focus our attention on (51) and (52) to establish the relationship between kinetic and kinematic quantities and ultimately establish the identity of the tensor . Various attempts have been made to develop constitutive models for the peridynamic theory [112,[114][115][116][117][118][119][120]. These material models can broadly be grouped into bond-based, ordinary state-based and non-ordinary state-based models. Let W _ ∶ V → ℝ be the peridynamic strain energy density, then generally, for an elastic material, the force density state can be expressed [112] [112] will be used in this study. This material model is based on the non-ordinary state-based framework that allows for the adaptation of the classical continuum model into peridynamics. A key motivation for this choice is to take advantage of decades of development and calibration of the classical model.
In the correspondence model, W _ is assumed to be equal to the strain energy density Ω( ) ∶ L 2 → ℝ from the classical theory where the deformation gradient from the local theory is approximated by its nonlocal counterpart given by (55). So that (62) becomes Evaluating the Fréchet derivative in (63) is shown [112] to result in the expression where P is the first Piola stress tensor. Comparing (64) to (42), it can be deduced that the force state _ is a state expanded from the tensor Comparing (51) and (65) shows that Thus, in compact notation, the state-based peridynamic balance of linear momentum can be written as To complete the definition of the nonlocal problem, we need to apply appropriate constraints to certain region of the problem domain. To this end, let the interaction domain Ω I = Ω c be the constrained volume. Let Ω c splits into two disjoint subdomains Ω cd and Ω cn such that Ω cd ∩ Ω cn = ∅ and either of Ω cd and Ω cn could be an empty set. Ω cd is the subdomain where Dirichlet boundary condition is applied and Ω cn is the subdomain where Neumann boundary condition is applied. Analogous to the boundary value problem of the classical local theory, constraint on the solution u of (49) over Ω is applied as follows: a given function value g d is prescribed on the solution over Ω cd such that To prescribe the Neumann-type constraint, recall that in the classical boundary value problem, this involves prescribing a traction or flux density • n over the traction boundary. From (18) and the discussion that follows, the nonlocal flux density over Ω In is given by Let g n be a given function value of the flux density over Ω cn , the Neumann constraint can be stated as The presence of the second-order time derivative of the solution u in (67) means in addition to the boundary constraints (68), (69), initial conditions also need to be specified. The initial condition involves prescribing the initial values of the solution and its first derivative. Let u I and u I be the initial values of u(x) and u(x), respectively, then and are the initial conditions. So that (67)-(71) gives the complete definition of the nonlocal problem: Notice that the constraints (68) and (69) are prescribed over domains Ω cd and Ω cn that have positive volume in ℝ n . This contrasts with the classical local theory where constraints are applied on domains that have zero volume. For this reason, (72) is referred to as initial volume constraint problem.

Discretization of the Peridynamic model
As can be seen from (72), the governing equation of motion in Peridynamics gives rise to a continuum model. To be amenable to computer implementation, different numerical approximation schemes have been proposed such as the meshfree method [121,122], the collocation methods [123,124] and methods based on finite-element mesh [125,126]. Due to its simple implementation algorithm and relatively low computational cost, the meshfree method suggested in [121] is the most widely used [127] and is the preferred method in this work for these same reasons. In this approximation method, the discrete form of (52) is and N is the number of nodes in the neighbourhood of node i.

Peridynamic homogenization theory
The homogenization procedure developed in this work is a two-scale scheme: a microscopic scale represented by an RVE, and a macroscopic scale represented by a homogeneous effective medium. The constitutive law of the microscale model is assumed to be explicitly known at every point of the micro-domain while the constitutive law of the micromodel is not known everywhere. The objective is to retrieve a constitutive law of the macroscale substitute medium from a numerical solution of an initial volume constraint problem (IVCP) at the level of the underlying microstructure. In this multiscale framework, an RVE is assigned to each integration point of the macro-continuum. A peridynamic equilibrium solution of the RVE is sought using boundary condition generated by the macroscale deformation gradient. The solution of the RVE IVCP yields the microscale stress field which is then homogenised to produce macroscale stresses and associated material tangent tensor. The coupling of the micro-and macroscale is achieved through averaging relationships and the energetic equivalence statement of the Hill-Mandel micro-homogeneity condition.

Effective material constants
Consider a heterogeneous medium o with characteristic size of heterogeneities to be l hetro . Momentarily, let this medium be replaced by a homogeneous 'effective' medium h . The original heterogeneous medium is the microscale, and the geometrical arrangement and material characteristics of the heterogeneities constitute the microstructure while the effective medium is the macroscale. Define a grid on h and let each point x on this grid be associated with a neighbourhood Ω s . Let Ω s be bounded by a region Ω c of positive volume inℝ n . Let a sample of h occupying the regions Ω s and Ω c be denoted as Ω 0 s and Ω o c , respectively. In addition, let a sample of o occupying Ω s be denoted as Ω o s as shown in Fig. 2. Now, define a grid on Ω s and let the position of points on this grid in the reference configuration be denoted as x . The grid associated with x is called the macroscale and the grid associated with x is the microscale. Let the characteristic lengths associated with the macroscale and the microscale be l macro and l micro , respectively. The morphology and material properties of the constituents of o in the microscale are called the microstructure of h . If  principle requires that the scale of the microstructure (or fluctuation of micro-field such as stress and strain) should be much smaller than the size of the RVE considered which in turn should be much smaller than the characteristic length scale of the macro-domain (or fluctuation of macro-field variables).
Let h be subjected to an affine deformation at its boundary. This will produce a homogeneous strain (for small deformation). This homogeneous strain will in turn generate a homogeneous stress field everywhere in h . For simplicity, we will assume linear elastic material response in both scales, then the material model that relates and given by the generalised Hooke's law: is the effective or homogenised constitutive law, where ℂ * and * are the effective stiffness and compliance tensors, respectively. At the microscale, the constitutive relation in each phase of the microscale is given by If the condition for the existence of the RVE is satisfied (henceforth, this condition will be assumed to be satisfied), then the microscopic deformation will be assumed to admit the following decomposition: where u * (x) is the displacement fluctuation due to the presence of the microstructure and * = G x u * (x).

Average theorems
The transition of mechanical properties from the microscale to the macroscale is achieved using volume average relations. Let be a quantity defined over a domain Ω . We denote the volume average of over Ω as  where V Ω is the volume associated with Ω.

Theorem 4.1 Nonlocal average stress theorem. Let be a heterogeneous body occupying the region Ω = Ω s
⋃ Ω c where Ω s is the region where solution is sought and Ω c is the boundary volume. We denote the average stress and average strain over Ω s as ⟨ ⟩ and ⟨ ⟩ , respectively. Let be in a state of static equilibrium when a constant stress tensor is applied on the boundary volume Ω c , then the volume average of the stress field in Ω s is equal to , that is Proof From (67), static equilibrium of the RVE in the absence of body forces requires the divergence of the Cauchy stress tensor in the case of small deformation to vanish, that is The Cauchy stress field in Ω s can be written as Taking the volume average of (81) and utilising (20) yields Considering the relationship D = −G * and utilising (80), then (82) reduces to       (78), the volume average of the strain field over Ω is given by This implies that the volume average of the strain field is completely defined in terms of the strain at the boundary volume.

Macrohomogeneity condition
For the averaged fields ⟨ ⟩ and ⟨ ⟩ to be admissible variables in the macroscale constitutive relation, the so-called Hill-Mandel macrohomogeneity condition [128] must be satisfied. The macrohomogeneity condition provides the basis for the substitution of an initially heterogeneous medium with a homogeneous one. This is achieved by prescribing an energetic equivalence between the heterogenous medium and the homogeneous substitution medium. Let the strain energy density of the underlying classical continuum material be then invoking the principle of constitutive correspondence allows us to write the strain energy of the peridynamic model as (85). The macrohomogeneity condition is stated as In other words, the condition requires that the average strain energy of the heterogeneous medium be equal to the strain energy density of the homogeneous medium. The condition under which (86) is satisfied for a peridynamic continuum material under constitutive correspondence will be established through a nonlocal analogue of the Hill's lemma. Proof We can write the scalar product of the average of the stress and strain tensor as Thus, we can rewrite (86) as which proves (87). It is obvious from (87) that for the Hill-Mandel condition, (86) to be satisfied will require that Thus, the satisfaction of the Hill-Mandel macrohomogeneity condition requires the integral in (89) to vanish.

RVE boundary volume constraints
In the classical continuum framework, the satisfaction of (89) is traditionally achieved in one of the following ways: (1) Voigt (or Taylor) assumption, (2) Reuss (or Sachs) assumption, (3) Homogeneous displacement, (4) Prescribed periodicity in displacement, (5) Homogeneous stress, and (6) Prescribed periodicity in traction. Methods 1-3 are categorised as deformation-driven approaches while methods 4-6 are categorised as stress-driven approaches. The task now is to establish the boundary requirements that will make the lemma (87) satisfy the Hill-Mandel condition in the nonlocal framework using the methods 1-6.

Voigt (Taylor) model
In this method [129], (89) is satisfied by assuming a homogeneous deformation of the form u i = x j ij in Ω . This implies a constant strain field (x) = in the RVE. Inserting this assumption into (75) yields (88) From (90), the implication of the Taylor (Reuss) assumption is that the homogenised or effective stiffness tensor is simply the volume average of the stiffness tensor of the constituents. It will also be noticed that utilising the Taylor assumption means, we can obtain the effective material properties without the need to solve the microscale peridynamic (RVE) problem.

Reuss model
In this model [130], (89) is verified by assuming a constant stress (x) = in Ω . If this assumption is inserted into (75) yields thus meaning that the effective compliance tensor is simply the volume average of the compliance tensor of the constituents. As with the Voigt model, the utilisation of the Reuss assumption allows the determination of the effective properties without recourse to solving the microscopic peridynamic (RVE) model.

Constant traction boundary volume constraint (CTVBC)
One way of satisfying (89) is to prescribe appropriate traction on the boundary volume Ω c . A traditional way of achieving this in the classical continuum framework is by applying the so-called constant traction boundary condition. In the nonlocal framework, this is achieved by imposing on the boundary volume Ω c , a constant traction generated by constant stress field Substituting (92) into (89) will vanish the boundary volume integral and, therefore, satisfying the Hill-Mandel condition (86).

Linear displacement boundary volume constraint (LDBVC)
This boundary condition is also referred to as homogeneous boundary condition in the literature. This boundary condition is achieved by applying appropriate displacement field to the boundary of the RVE that will varnish the gradient of the displacement terms of the integrand of the boundary volume integral (89). A traditional way of achieving this is to apply a linear displacement of the form Inserting (93) into (89) yields Thus, proving (93) satisfies the Hill-Mandel condition (86).

Periodic boundary volume constraint (PBVC)
This model is appropriate to model materials with periodic microstructure. The reference configuration of the RVE is assumed to be a geometric shape with even number of sides or faces for two-and three-dimensional problems, respectively. A square RVE is shown in Fig. 2 for twodimensional problems, with each pair i of the RVE boundary region assumed to be equally sized subsets, that is, there should be a one-to-one correspondence between points in Ω + i and Ω − i . In this method, a displacement field of the form (77) is applied on the boundary region such that for each pair of boundary points ( The difference in displacement between two corresponding boundary points x + , x − is then given by To achieve static equilibrium of the RVE, an anti-periodic stress field is applied in the boundary domain such that for each pair of points in Ω + i and Ω − i . Utilising (95) and (96), the Hill-Mandel condition is satisfied as follows:

Bounds for effective properties
Predictions from the Voigt and Reuss assumptions were shown in [131] to provide the upper and lower bounds for the effective material properties. That is, where ℂ * R , ℂ * V are the effective stiffness tensors due to Reuss and Voigt assumptions, respectively, and ℂ * is the effective stiffness due to any other method. As emphasised in [132], the corresponding entries of ℂ * R , ℂ * and ℂ * V do not necessarily satisfy (97), however, the corresponding diagonal terms and eigenvalues are shown to satisfy (97). It can also be shown that the inequality (97) hold true in terms of elastic constants such as bulk and shear moduli. Consider a given strain tensor kl at a point, the stress ij acting at the point for a linear 2D isotropic media undergoing small deformation is given by (96) and the fourth-order isotropic tensor ℂ ijkl can be written as where and are the bulk and shear moduli, respectively. If we define two fourth-order isotropic tensors I 1 and I 2 as then (98) can be written as It can be shown from (100) that the compliance tensor can be written as For the sake of brevity, we will adopt a symbolic notation that will allow us write (100) and (101) as Given a two-phase composite having constituent material constants 0 , 1 , 0 and 1 and volume fractions c 0 and c 1 such that 1 > 0 and 1 > 0 , then the volume average of the stiffness and compliance tensor can be written as Utilising (102) and (103) in (97), it can easily be shown that

It follows from (104) that
and where the left-hand sides of (105) and (106) give the Reuss lower bound while the right-hand sides give the Voigt upper bound, and * and * and the 2D effective bulk and shear moduli. It is noted that the distance between the Reuss lower bound and the Voigt upper bound is large and often does not give much information, a tighter bound is achieved using the Hashin-Shtrikman bounds [133]: and

Computational implementation of the PDCHT
To obtain the numerical solution of the RVE in the PDCHT framework, the RVE will be discretised following the procedure outlined in Sect. 3.3. Being a nonlocal problem, the RVE is subjected to appropriate volume constraints. In the numerical validation that follows in Sect. 5, the RVEs will be subjected to LDBVC and PBVC. Although computational algorithm to implement these . Fig. 3 Example square RVE showing corresponding boundary regions boundary conditions in the framework of finite-element analysis is well established and discussed by many authors [134][135][136][137], however, implementing them in a nonlocal boundary value constraint problem such as the RVE in the PDCHT framework require special treatment. This is particularly the case with the PBVC. To implement the PBVC in the PDCHT framework, the displacement-driven approach to homogenization is utilised in this work and thus, to determine the effective elasticity tensor, the first of (75) is employed. Combining this with (79) and (84) allows us to write the expression for the effective elasticity tensor as where the stress field ij in Ω s is obtained using discretised peridynamic equation of motion (73) and kl is the prescribed strain tensor on the boundary volume. For a twodimensional problem, ℂ * ijkl has six components. However, owing to its symmetric property, only three components are independent. Thus, to determine the components of ℂ * ijkl for a two-dimensional problem requires the application of three loading conditions that results in deformation modes which render all but one of the three independent components of the strain tensor to zero. For the purpose of this implementation, these deformation modes are given as where c is the magnitude of the prescribed strain tensor. These strains are then used to generate displacement in the boundary volume depending on the boundary condition used. In the case of LDBC, the displacement u generated at every node x i in the boundary volume after discretization of the RVE follows from (93) as . pair Ω − 2 and Ω + 2 , respectively, it will be noticed that each corner volume is shared by two facial volumes and applying (102) to nodes within these boundary volumes will lead to over constrained boundary. To eliminate this problem, the following translational periodicity is imposed at the corner boundary volumes: if we write

Fig. 4 RVE geometry showing various configuration
and Ω + 4 = Ω c 4 , then volumes Ω c 1 and Ω c 2 are assumed to be images of Ω c 3 under horizontal and vertical translational symmetry, respectively, while Ω c 4 is the image of Ω c 3 under combined horizontal and vertical symmetry so that the following relative displacement constraints are imposed between nodes in the corner boundary volumes: To eliminate translational rigid body motion of the RVE, nodes within the corner volume Ω c 3 are constrained. The determination of the effective elastic tensor in both problems proceeds under the assumption of small deformation and plane stress. This allows the use of infinitesimal strain tensor and Cauchy stress tensor directly in (64). Once the Cauchy stress field in the RVE is obtained, (98) is used to recover the effective elasticity tensor.

Validation of the homogenization scheme
Having established and justified the peridynamic correspondence homogenization theory (PDCHT) in Sect. 4, numerical examples are presented in this section to benchmark the scheme against the Reus-Voigt and Hashin-Shtrikman bounds discussed in Sect. 4.2.3. Prediction from the PDCHT would also be compared against predictions from (116)  other classical mean-field homogenization methods such as the Eshelby dilute estimate [138] and the Mori-Tanaka method [139]. In addition, the result from the PDCHT will be compared to result obtained from computational homogenization based on the Finite-Element Analysis. After validation of the PDCHT strategy, the framework will be used to predict the effective properties given elliptical inclusion to observe the influence of inclusion shape on the effective properties predicted by the method. All materials considered throughout this section are assumed to be two-phased consisting of a matrix phase and a stiffer fibre phase (inclusion). Properties associated with matrix will be denoted with the superscript (m) while those associated with the fibre phase will be denoted with superscript (f ) . Both matrix and fibre phases are assumed to be isotropic under isothermal linear elasticity.
To pursue the objectives of validating the proposed method, three numerical examples will be considered. Figure 4 shows the RVE configurations to be considered in this section, and the properties of the constituent materials for the RVEs to be considered in the numerical examples are given in Table 1.

Comparing the PDCHT results with bounding theorems and other established methods
In this example, effective properties predicted from the PDCHT will be compared against computational result from the bounding theorems of Reuss, Voigt and Hashin-Shtrikman, the mean-field methods of Eshelby and Mori-Tanaka as well as the full-field method based on the finite-element analysis as done by the authors. The material is assumed to be a glass in epoxy-matrix composite with properties given in Table 1 under a plane strain condition. The RVE geometry is assumed to consist of epoxy matrix with a circular glass fibre centrally placed as shown in Fig. 4a. The problem is solved over the range of all admissible fibre volume fractions 0-100%. Solutions will be sought considering LDBVC and PBVC. Results from the computations are presented as follows. The predicted evolution of the effective elastic stiffness tensor and corresponding effective elastic constants from PDCHT, the bounding theorems and other established methods are presented in Figs. 5 and 6, respectively, for the case of LDBVC. Similar analysis with the same RVE subjected to PBVC is conducted and the results are presented in Figs. 7 and 8 representing the evolution of the effective stiffness tensor and elastic constants, respectively.
Prediction from PDCHT of the effective stiffness tensors presented in Figs. 5 and 7 lie within the Reuss-Voigt bound as well as the tighter Hashin-Shtrikman bound thus satisfying (97). Similar agreement with the bounding theorems is observed with prediction of the effective bulk and shear moduli as presented in Figs. 6 and 8 thus satisfying (105) and (106) for the Reuss-Voigt bounds and (107)-(108) for the Hashin-Shtrikman bounds. This is true for both the predictions under LDBVC and PBC. Another consequence of the bounding Eqs. (97), and (105)- (108) is that the elastic modulus from any proposed homogenization theory is predicted to lie within the Reuss-Voigt and Hashin-Shtrikman bound. The evolution of the effective elastic modulus obtained from the proposed theory indeed lies within these bounds as shown in Figs. 6c and 8c for LDBVC and PBVC, respectively.
The result of prediction from the PDCHT is also compared to predictions from the mean-field homogenization methods of Eshelby and Mori-Tanaka as well as a fullfield computational method based on FEM solution of the RVE. The results from these methods are also presented in Figs. 5, 6, 7 and 8. The predictions from both the PDCHT and FEM comply with the bounding theorems and compared to predictions from the Mori-Tanaka method, the PDCHT prediction gives an upper estimate of the effective properties. It is worthy to note that the prediction in this example from the Mori-Tanaka method coincide with those from the Hashin-Shtrikman lower bound. This is the case in some situations and has been reported in the literature [González]. Predictions from the Eshelby dilute method agrees with results from the PDCHT only for very low fibre volume fraction. This is an expected trend as the dilute method is expected to give reasonable predictions only for very low (dilute) fibre volume fraction.
Since the fibre is assumed to be of circular cross-section, it can be shown that the maximum volume fraction that can be achieved with a perfectly circular fibre crosssection is 78.54% . Beyond this volume fraction, the circular geometry of the fibre cross-section degenerates and thereby causes a change in the morphology of the RVE. This change in morphology is reflected in the results predicted from the PDCHT and FEM solutions of the RVE and expectedly not the bounding theorems as shown in Figs. 5, 6, 7 and 8. A comparison of the predicted effective elastic modulus obtained using both LDBVC and PBVC as shown in Fig. 9 shows that the prediction using the LDBVC provides an upper estimate of the two boundary conditions at least within the small deformation regime. This is a tested and proven result from the literature [132,140,141].

Comparing the PDCHT results with results from published works
The objective in this example is to compare predictions from the proposed PDCHT with experimental results from [142]      and computational predictions from the following references: [105,108,[143][144][145]. The RVE is assumed to be a square array of circular shaped boron fibre placed in the centre of aluminium matrix as shown in Fig. 4a. The properties of boron and aluminium are given in Table 1. The RVE volume-constrained problem is solved under the assumption of plane stress and PBVC. Prediction of effective elastic properties of the composite system from the PDCHT is presented in Table 2 alongside results from some published references. Analysis of the results shows that the prediction from PDCHT provides the closest correlation to the experimental result from [142] in the estimate of the effective elastic and shear moduli. Compared to other computational methods, prediction of these moduli from the PDCHT gives the lower estimates. However, the effective Poisson's ratio 12 = 0.185 predicted by the PDCHT is markedly different from the experimental result but agree well with predictions from the FEM, OSB-PDHT, FEM, PD UC and LTEHOT.

Effective properties of RVE with elliptical fibre inclusion
In this example, two RVEs, the first with circular and the second with elliptical fibre inclusions as shown in Fig. 4a and b, respectively, will be considered. Two composites will be considered. The first is a glass in aluminium-matrix composite and the second is a graphite in aluminium-matrix composite with material properties as given in Table 1. The stiffness ratio (or phase contrast) ( ∶= E (f ) ∕E (m) ) of the first material is 1.07 while that of the second material is 3.44 . The objective in this example is to briefly demonstrate the capability of the PDCHT in capturing the effect of inclusion shape and phase contrast on the effective behaviour of materials using the LDBVC. The evolution of the normalised effective elasticity modulus (effective stiffness ratio) * i = E * i ∕E i(m) with respect to fibre volume fraction c (f ) for RVEs with both circular and elliptical inclusion for the two materials alongside results from FEM simulation of the same problems are presented in Fig. 10. In Tables 3, 4, 5 and 6, the evolution of the elastic modulus in directions 1 and 2 as well as the percentage difference between predictions from PDCHT and FEM are presented.
From the result presented for predictions from both PDCHT and FEM, the effective material behaviour as represented by the evolution of the effective elastic modulus for the RVE with circular inclusion shows linear behaviour for low material phase contrast ( = 1.07) and weak nonlinear behaviour for higher material phase contrast ( = 3.44) . In the case of RVE with elliptical inclusion, the result for lower material phase contrast shows linear behaviour while the effective material behaviour at higher material phase contrast shows strong nonlinear behaviour. The amplified nonlinearity is because the elliptical inclusion introduces an anisotropy in the microgeometry of the composite. In addition, the prediction of effective material properties for RVEs with circular inclusion from Tables 2  and 3 shows an isotropic effective material response while predictions for RVEs with elliptical inclusion shows an orthotropic effective material behaviour. Comparing the results from the PDCHT and FEM, it is noted that at lower material phase contrast, the predicted effective behaviour from the two methods shows good agreement for both inclusion shapes especially at low material phase contrast. The maximum percentage difference in the estimated results from these methods is 0.421% for circular inclusion and 0.45% for elliptical inclusion. However, at higher material phase contrast ( = 3.44) , the maximum percentage difference raises to 6.55% for circular inclusion and 9.2% in the case of elliptical inclusion, with the predictions from the PDCHT yielding more nonlinear behaviour with increasing material phase contrast. This nonlinear effective material behaviour is an expected prediction and has been reported in literature [146]. In addition, the difference in result between the PDCHT and FEM predictions at high material phase contrast also did not come as a surprise because it has been reported in the literature [147,148] that marked difference between different homogenization approaches at high stiffness ratio have been observed.

Conclusion
This paper presented a first-order homogenization theory in the framework of the non-ordinary state-based correspondence model. The development of the theory is set on a rigorous mathematical framework consistent with the nonlocal nature of the peridynamic theory. Using this homogenization theory, it is now possible to obtain microstructure informed properties of materials for use at the macroscale within the framework of peridynamic modelling. The proposed nonlocal homogenization theory is validated by solving benchmark problems and comparing the results with those obtained using the Reuss, Voigt and Hashin-Shtrikman bounding theorems, the mean-field methods of Eshelby and Mori-Tanaka as well as the finite-element method. These results are shown to comply with the bounding theorems as well as agree with results from the mean-field and full-field methods mentioned above.
Comparison of predictions from the proposed homogenization framework with results from the literature shows good agreement. Lastly, the PDCHT has been shown to be capable of capturing interesting material behaviour some of which have been reported in the literature such as the effect of change in the effective material property when the morphology of the RVE changes as illustrated in the first example. This effect was demonstrated more elaborately in the third example by comparing the prediction for RVEs with circular and elliptical shaped inclusion. The proposed method was also shown to be able to capture the effect of material phase contrast on the effective behaviour of the composite system.
The advantage of PDCHT theory over homogenization frameworks based on the classical continuum theory derives from the strengths of the peridynamic theory. This framework is, therefore, especially useful in circumstances involving evolution of the microstructure or problems in which nonlocal interaction plays important role in the overall response of the heterogenous media. Another advantage that can be leveraged with this development is that because the peridynamic correspondence model uses familiar quantities from the classical continuum theory, Once the effective material tangent is obtained, we are free to use either of peridynamic theory or the classical continuum theory to characterise the macroscopic response of the medium. Where the peridynamic theory is utilised at the macroscale, this result in a standard multigrid method we will call the PD 2 method. In the case where the classical theory is utilised, this result in what is referred to in the literature as Heterogenous Multiscale Method (HMM) [149]. In this case, numerical schemes such as the finite-element method or the finite difference method can be utilised to solve the macro-model.
Since one of the greatest strengths of the peridynamic theory lies in its capability in handling strong discontinuity in the response field of a system, this homogenization theory will be applied to problems with evolving microstructure such as micro-crack coalescing and propagation in a future work.
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:// creat iveco mmons. org/ licen ses/ by/4. 0/.