Peridynamic computational homogenization theory for materials with evolving microstructure and damage

This study aims to establish a framework for multiscale assessment of damage for materials with evolving microstructure based on a recently proposed peridynamic computational homogenization theory. The framework starts with replacing a material with complex microstructure with a constitutively equivalent material that is microstructurally homogenous. Constitutive equivalence between the original and the substitute materials is achieved through enforcing strain energy equivalence via the so-called nonlocal Hill’s lemma. The damage law is obtained by numerically solving boundary volume constraint problem of an RVE. The result from the analysis of the RVE problem was compared with the previously published result to establish the validity of the proposed framework. The comparison shows good agreement between result obtained using the proposed framework and those reported in the literature.


Introduction
Peridynamics (PD) [1] has been used to solve a wide range of engineering problems  since its introduction, and particularly, it is increasingly used to study the behaviour of composite materials [23][24][25][26][27][28][29][30][31][32]. When modelling composite materials with a complex microstructure resulting from manufacturing processes and oftentimes random distribution of constituent particles [33], it is a common practise to replace highly heterogeneous properties with effective material properties to avoid the costly computational requirement of numerical simulation that will directly resolve information from all relevant scales for the entire computational domain. Apart from the effective properties of nonhomogeneous materials like composites, there is also an increasing need to incorporate microscale damage processes into phenomenological models of materials that would otherwise be deemed homogeneous. Again, a straightforward numerical simulation that captures all important microscale mechanisms will be computationally expensive.
These considerations necessitated the development of different classes of multiscale modelling framework which were designed to allow for the simultaneous use of many models at various scales of resolution to characterize the behaviour of a material. Conceptually, depending on the type of problem, these multiscale frameworks can be grouped into three broad categories, viz.
1. Type 1. Interactions in a subscale region is lumped into a point of the macrodomain. 2. Type 2. The computational domain is divided into subdomains in which either.
a. Coarse mesh representing the macro-continuum and fine mesh representing material at the subscale is used in different subdomains, or b. Different models are used in different subdomains to resolve relevant information.
3. Type 3. The macroscale domain is discretized, and every macro-point is embedded with a subscale material element.
Coarse graining or model-order reduction are terms used to describe the Type 1 multiscale technique. The response of a discretized macro-continuum is tracked using a reduced number of points on the macro-grid in this framework. In 1 3 other words, this framework enables the solution of a highfidelity model with fewer degrees of freedom. Two multiscale strategies from this category have been proposed. Silling [34] proposed a coarsening method in which a series of successively coarsened body is derived from a more detailed material description. The coarsening method was extended for two-dimensional implementation in [35]. Another coarse graining methodology for PD was proposed in [36] based on static condensation algorithm. The order of the model is reduced by ignoring the inertial component of the equilibrium equation and expressing the discarded degrees of freedom in terms of the retained degrees of freedom. This method was employed in [37] to coarsen PD heat transport model. See [38] for a comprehensive treatment of type 1 methods.
In the type 2 multiscale methodologies, the same or different models focusing on different physics or scales of resolution is/are applied to finite subregions of the computational domain. Multiscale methodologies falling under this category are called concurrent multiscale methods. Two broad categories of concurrent methodologies have been proposed for PD, namely monomodel and multimodel methods. PD is used to represent the whole computational domain in monomodel approaches [39][40][41][42][43][44], and multiscale capability is achieved by appropriately refining the grid to resolve features at multiple length scales. Multimodel concurrent methods are hybridizing frameworks that aim at creating coupling strategies that permit peridynamics to be implemented alongside other modelling frameworks such as classical continuum or atomistic modelling theories. The goal is to take advantage of the benefits that each of the connected models has to offer. Multibody concurrent techniques can be characterized as kinematic-based [45][46][47][48], force-based [49][50][51][52][53][54][55], or energy-based [56][57][58][59][60][61][62] methods depending on the coupling mechanism utilized.
Methodologies that fall under the type 3 multiscale framework are generally referred to as the hierarchical multiscale method in the literature. The objective of material characterization using this multiscale framework is to evaluate the effective material properties at the macroscale based on mechanisms of deformation and damage occurring at a subscale. (Henceforth, the subscale will be referred to as the microscale.) The key assumptions made in developing these methods are that while inherently heterogeneous at the subscale, the macro-continuum is assumed to be effectively homogeneous. This assumption is deemed justified where the principle of scale separation is achieved. The process by which a heterogeneous material is replaced by a constitutively equivalent homogeneous material is called homogenization; hence, methods in this category are also called homogenization methods. The assumption of macrohomogeneity is usually justified either in a statistical sense or based on periodicity of the subscale structure. Where statistical justification is used, a representative volume element (RVE) is employed to represent the microscale material. The RVE is the smallest sub-volume of a statistically homogeneous microstructure that has the same phase volume fractions and statistical distributions of constituents and other microstructure characteristics as the rest of the material and responds to homogeneous displacement or traction boundary conditions in the same way as the rest of the assemblage. On the other hand, if macro-homogeneity is justified through periodicity such that subscale structure can be represented by a stack of repeating cells, then a single of such cells is called a repeating unit cell (RUC). It is worth noting that because a periodic structure is a particular case of random structure, methods based on the notion of RVE can be utilized to homogenize a medium with periodic microstructure.
A few homogenization techniques have thus far been proposed and implemented within the PD modelling framework. A method of cell was proposed in [63] in which micro-macro coupling is achieved through a strain concentration tensor which relates the macroscale deformation to the microscale deformation. The strain concentration tensor is a complex function of the geometry and constitutive properties of the constituent materials of the unit cell and is defined to express the fluctuating micro-strain in terms of the external average strain. The components of this tensor are obtained by solving an equilibrium problem of the unit cell subjected to periodic boundary conditions. The method of cell relies fundamentally on the assumption that the heterogenous medium has a periodic microstructure.
A family of PD homogenization methods based on strain energy equivalence have been proposed and utilized in [64][65][66][67][68]. The key idea in these methods is that two materials with different microstructure are taken to be constitutively equivalent if the strain energy stored in both materials under affine deformation is the same. While the method proposed in [67,68] was proposed for application within the BBPD modelling framework, the method proposed in [64] was developed for application in the framework of the ordinary state-based peridynamics and is called ordinary state-based peridynamic homogenization theory (OSBPDHT). It was utilized to study the effective behaviour of elastic materials with evolving microstructure.
In [65], a mathematically rigorous homogenization framework that is compatible with the nonlocal nature of PD was proposed. This method is called peridynamic computational homogenization theory (PDCHT). The PDCHT is a two-scale micro-macro homogenization technique in which the macroscale constitutive relation is generated from explicit solution of a microscale nonlocal volume constraint problem. Nonlocal analogues of the stress and strain average theorems, as well as the Hill-Mandel macrohomogeneity condition, were derived to justify the coupling between the two scales. To establish the validity of the method, solution to numerical problems was sought and compared with the previously published results.
Motivated by the development of the PDCHT in [65], the contribution of this work is to extend the application of the PDCHT to a class of materials with evolving microstructure. In other words, this study is interested in examining the influence of evolving microstructure on the effective response of materials at the macroscale. The effective behaviour represents the mechanical state of a point at the macroscale and is defined to be the average behaviour of the microstructure over a finite micro-subregion. The material properties arising from this state then give the effective material properties at that point. Because the microstructure informs the mechanical condition of the macro-continuum point, an evolving microstructure will thus result in changes in material characteristics at the macroscale. Evolution of the microstructure refers to energy dissipating processes occurring at the microscale such as nucleation, propagation, and coalescing of crack. The irreversible changes brought about in a material by these microscale processes that result in degradation of the material is referred to as damage.
To proceed with the aim of this contribution, we present in Sect. 2 the formulation of the homogenization framework in the context of evolving microstructure. Section 3 presents detail procedure for computational implementation. Benchmark problems are solved in Sect. 4 and the results are compared with results from [69] and the OSBPDHT proposed in [64]. It is worth noting that the difference between the method proposed in this contribution and the OSBP-DHT is that OSBPDHT is developed in the framework of the OSBPD while the present contribution is based on the NOS-BPD framework. The benefit we aim to gain with our proposed framework derives from the strength of the NOSBPD, which allows for a more generalized material model than the OSBPD permits. Another advantage of the NOSBPD that is exploited in this work is that it uses the familiar stress and strain concepts from classical continuum theory, therefore making hierarchical coupling of PD and CCM models across scales a more seamless endeavour.

Peridynamic computational homogenization
Consider an elastic peridynamic body Ω with a heterogeneous microstructure subjected to prescribed displacement on some or all its external boundary volume Ω Γ as shown in Fig. 1. Let the microstructure be characterized by microconstituents and the presence of microvoids and/or microcracks. Let the micro-constituents be assumed to possess homogeneous properties, and let Ω satisfy equilibrium, then the non-ordinary state-based peridynamic boundary volume constraint problem is given by the following set of equations:

Fig. 1 Macro-and micro-grid
Equation (1) is the NOSBPD correspondence model. Given a bounded open domain Ω ∈ ℝ n , under the state-based PD framework, A continuum point x ∈ Ω interacts with an infinite number of other points inside its influence domain. If this domain is assumed to be a ball B (x) of radius > 0 centered at x , then is called the horizon of x , such that: where the set of all x � ∈ B (x) is called the family of x . The relationship between two interacting material points x and x ′ is called a bond and the relative position vector | | = |x ′ − x| in the undeformed configuration is called the bond length. (1) is referred to as the force state. The force state is the force vector per unit volume squared between x and x ′ . Generally, states are functions specified on bonds in B (x) . Another vector state that is relevant in the state-based PD formalism is the deformation state. The deformation state is a function acting on each bond length = x ′ − x for every x ′ in B (x) . The value of is the image of the bond in the deformed configuration given as: 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 . In (1), the tensor is given by To specify the domain H of a state, let > 0 be the horizon of a point x in a PD body , then: defines the set of bonds that exist between the point x and all x � ∈ B (x) . The nonlocal infinitesimal strain tensor in (1) is given by where in (6) is the nonlocal deformation gradient given by The variable u in (7) represents the displacement and is a weight function.
To numerically solve (1), the domain Ω must be discretized into a grid and a discrete form of (1) applied to each point on the grid. Different discretization strategies have been suggested for the PD equation of motion. These methods include the meshfree method [70,71], collocation methods [72,73], and methods based on finite element mesh [74,75]. The meshfree technique proposed in [71] is the most extensively used and is the preferred method in this study for reason of its easy implementation and relatively cheap computing cost. In this approximation method, the discrete form of (1) is and N is the number of nodes on the grid in the neighbourhood of node i.
For many practical problems, discretizing the whole macrodomain with a grid spacing that is capable of directly resolving relevant microscale characteristics is at best not computationally efficient and sometimes simply not possible given the present limited availability of computational resources. For this reason, a multiscale modelling framework based on PDCHT is proposed in this study to account for the effect of microstructure on the macroscale response without having to explicitly resolve the detail of the microstructure. This is accomplished by substituting the original body Ω having a heterogeneous microstructure with a constitutively equivalent material Ω * with yet to be known "effective" properties. Determining the effective material properties of the substitute homogeneous material Ω * will be achieved by considering RVEs from both Ω and Ω * . Let the RVE from Ω be denoted as Ω RVE and the RVE from Ω * be denoted as Ω * RVE . The starting point for the energy-based family of homogenization theories upon which the proposed framework depends on is that constitutive equivalence between and Ω * is achieved if under the same macroscopic deformation, the same strain energy is stored in both Ω RVE and Ω * RVE . This is the statement of the celebrated Hill-Mandel's macrohomogeneity condition which for PD correspondence model as is with the classical continuum theory is given as If energetic equivalence is achieved between Ω RVE and Ω * RVE as expressed in (9), then Ω * RVE can be replaced with Ω RVE in Ω * . If this is the case, then if an affine deformation is applied to Ω * at its boundary volume, and linear elastic behaviour is assumed, this will result in a homogenous strain * which will consequently generate a homogeneous stress field * everywhere in Ω * including Ω * RVE ⊂ Ω * . On the other hand, if we replace Ω * RVE with Ω RVE in Ω * , then the affine deformation will produce a nonhomogeneous strain field (x) which in turn will produce a nonhomogeneous stress field (x) in Ω RVE . In the NOSBPD correspondence material model the relationship between (x) and (x) and between * and * are, respectively, given by Hooke's law as and In this homogenization framework, it is assumed that the structure, constitutive law, and mechanical properties of the constituent materials in Ω RVE are explicitly known. This means ℂ(x) is known and with a well-posed boundary volume constrained problem (x) and (x) in (10) can be computed for every point x in Ω RVE . This permits for the reduction of the computational domain of (1) from Ω to Ω RVE . Once the (x) and (x) are obtained for every point x in Ω RVE , the task of determining the effective material tensor ℂ * is achieved by establishing relationships between stress and strain fields (x) and (x) in Ω RVE and their counterparts * and * in Ω * RVE . This relationship is established via the average stress and strain theorems.
Consider a heterogeneous body occupying the region Ω = Ω RVE ⋃ Ω RVE Γ , such that Ω RVE is the region where solution is sought and Ω RVE Γ is the boundary volume of Ω RVE . Let the average stress and average strain over Ω RVE be denoted as ⟨ ⟩ and ⟨ ⟩ , respectively. Let be in static equilibrium when a constant stress tensor is applied on Ω RVE Γ . The nonlocal average stress theorem [65] allows the following preposition to hold true: if the heterogeneous body Ω RVE is subjected to a traction boundary condition generated by the constant stress tensor , then irrespective of the complexity of stress field ( (x)) within Ω RVE , the stress averaged over Ω RVE is the same as . If is chosen to be * , then the statement of the average stress theorem can be written mathematically as Similarly, the nonlocal average strain theorem allows the following preposition to be true: if Ω RVE is subjected to applied displacement on its boundary volume Ω RVE Γ which is generated by a constant strain tensor such that u 0 = x ∀x ∈ Ω RVE Γ , then irrespective of the complexity of the infinitesimal strain field ( (x)) generated in Ω RVE , the infinitesimal strain field averaged over Ω RVE is the same as . Again if is chosen to be equal to * , then the average strain theory can be mathematically stated as What the average statements (12) and (13) allow us do is to couple the "micro-fields" (x) and (x) defined in the microscale with "macro-fields" * and * defined on the macro-continuum scale. The next very important task in developing this homogenization framework is to determine the admissible fields * and * that will satisfy the energy equivalence criteria (9). This task is achieved by determining the condition under which * and * will satisfy the so-called nonlocal Hill's lemma which was proved in [65] to be where S s x k is a nonlocal symmetric tensor operator (see [65] and [76] for explanation on nonlocal integral operators). It is obvious that for the Hill's lemma to satisfy the macrohomogeneity condition (9) will require the integral in (14) to vanish. In other words, a sufficient condition to satisfy the energy equivalence criteria is that There are many ways of satisfying (15). However, in the context of peridynamic computational homogenization framework [65], this has been achieved in one of three ways of prescribing boundary volume constraint to Ω RVE .

Boundary volume constraints on RVE
In the framework of PDCHT [65], (15) can be satisfied by either prescribing a constant traction boundary volume constraint (CTVBC), linear displacement boundary volume constraint (LDBVC), or periodic boundary volume constraint (PBVC).

Constant traction boundary volume constraint (CTVBC)
This boundary volume constraint is accomplished by imposing a constant traction on the boundary volume Ω RVE Γ created by a constant stress field of the form.
Substituting (16) into (15) will make the boundary volume integral evaluate to zero and therefore satisfying the Hill-Mandel condition (9).

Linear displacement boundary volume constraint (LDBVC)
This boundary volume constraint is accomplished by introducing a displacement field to the RVE's boundary volume that vanishes the gradient of the displacement terms in (15). One way of achieving this is to apply a linear displacement of the form: Inserting (17) into (15) yields Thus, proving that applying a LDBVC of the form (17) will guarantee that the Hill's lemma (14) satisfies the Hill-Mandel macrohomogeneity condition (9).

Periodic Boundary Volume Constraint (PBVC)
If the condition for spatial periodicity is taken to hold true in the microscale, then this means a mechanical field within the RVE can be split into its average and fluctuating components. In this respect, the displacement and strain fields within the domain Ω = Ω RVE ⋃ Ω RVE Γ can be written as

18) u(x) = * x +ũ(x) and (x) = * +̃
To apply the periodic boundary constraint on the RVE boundary volume, the boundary volume is divided into a positive part Ω + RVEΓ and a negative part Ω − RVEΓ and a displacement of the form (18) is applied such that for each pair of boundary points (x + ∈ Ω + RVEΓ , x − ∈ Ω − RVEΓ ) It is obvious from (18) and (19) that the difference in displacement between two corresponding points x + , x − is given as An anti-periodic stress field is applied in the boundary domain to attain static equilibrium of the RVE, so that for each pair of boundary points x + ∈ Ω + RVEΓ and x − ∈ Ω − RVEΓ , we have It is easily verified that using (19) and (21) in (15) satisfies the Hill-Mandel macrohomogeneity condition (9):

Peridynamic implementation
The strategy for implementing the homogenization framework developed in the preceding sections can be described by the following steps.
1. The originally heterogeneous peridynamic body is assumed to be macroscopically homogeneous and discretized using for example the meshless method [71] into nodes which when taken together forms a grid. The nodes and grid in this scale are denoted, respectively, as the macroscopic nodes and macroscopic grid. 2. An RVE(P) is assigned to each node P on the macroscopic grid. The RVE is chosen to adequately describe all relevant microstructural morphology of the original heterogeneous material. 3. For each macroscopic node P, compute the evolution of the local macroscopic strain * . This macroscopic strain is then used to impose appropriate constraint to the boundary volume of the assigned RVE. Since the microscale RVE boundary constraint depends only on the first-order gradient of the macroscopic displacement field, this method is also called first-order homogenization. 4. The microscale volume constraint problem is solved using (1) to yield deformation u(x) , strain (x) and consequently the stress (x) field within the RVE(P). 5. The stress field (x) is used in (12) to compute the average stress, * over the RVE(P). 6. For each macroscopic node P, substitute * (P) and * (P) into (11) to compute the macroscopic stiffness tensor ℂ * (P).
This algorithm is implemented into a first-order Peridynamic computational program in MATLAB 2021.

Numerical implementation of the first-order homogenization
We present two numerical problems involving crack growth and accumulation in this section to demonstrate the validity of the proposed homogenization framework. The effective elastic tensor is determined in both problems based on the assumptions of small deformation and plane stress condition, thus justifying the utilization of the infinitesimal strain tensor and the Cauchy stress tensor in (1). To accomplish this objective, the benchmark results from [64,69] are used. Since damage accumulation due to nucleation and propagation of cracks is known to induce anisotropy in the elastic properties of materials, we provide a brief remark on a suitable measure to quantify the crack induced anisotropy on the effective properties of the example material. Different measures [77][78][79] exist that seek to quantify elastic anisotropy of materials. However, the method proposed in [80] is preferred in this contribution due to its simplicity and most importantly its suitability for two-dimensional problems. In this method, the elastic anisotropy index A for a twodimensional crystal is given by where C ij and S ij are the components of the elastic stiffness and compliance tensors, respectively. In (22), A yields a value of zero for the limiting case of perfect isotropy while deviation from zero defines the degree of elastic anisotropy.

Material softening due to crack propagation
In this example, the proposed homogenization scheme will be utilized to study the phenomenon of material softening due to the presence and interaction of cracks at the microscale. To this end, two scenarios are investigated and the RVEs considered are shown in Fig. 2. Each RVE is taken to be a square of side length W = 1m and comprise of a matrix material with an elastic modulus E = 2 × 10 9 Pa and (22) a Poisson's ratio v = 0.3 . In the first scenario, the RVE in Fig. 2a is utilized, which consists of a single horizontally propagating crack whose origin is at the centre of the RVE. The RVE for the second scenario is given in Fig. 2b and consists of two coalescing cracks. In the first case, the length of the crack at any given instant of time is L , while in the second scenario, the length of each crack is L∕2 . In both scenarios, the initial length of the crack is taken to be L = 0 and the final total length is L = 1m.
We remark that although one of the key advantages offered by peridynamics is in the modelling of autonomous crack propagation, in this contribution, the propagation of the crack is controlled to be consistent with the propagation in the referenced contributions [64,69]. To achieve the desired crack length consistent with crack propagation in the referenced works, the analyses in this example are conducted under a quasi-static loading assumption through multiple load steps. The realization of each load step is represented by a specific length of crack.
The simulation is done over crack propagation processes spanning twenty-one discrete crack lengths in the interval [0 − 1] and [0 − 0.5] for the first and second scenarios, respectively, which translate to a crack growth of 0.05 m for each load step. The distance r between the origin of the two cracks in Fig. 2b is given as 0.25m , thus allowing the cracks to coalesce at L = 0.5m.
The effective tangent matrices of the material were evaluated as its microstructure evolved due to the single propagating crack and two coalescing cracks, respectively, following the technique outlined in Sect. 3. Evolution of the C 22 component of the effective tangent stiffness matrices for the case of a single propagating crack is presented as Fig. 3a and for the case of two coalescing cracks is presented as Fig. 3b. These results are compared with results from [64,69]. The results from the single propagating crack scenario indicate a smooth relationship between material softening and crack growth. However, in the case of the two cracks, an abrupt change in the relationship is observed as the two cracks coalesce. This phenomenon is attributed to interaction between the coalescing cracks. It was noted in [69] that when multiple crack coalesce, this leads to significant increase in a damage energy release rate which lead to amplification of damage which in this case is characterized by the sudden change in material softening as observed. In both cases, it is observed from Fig. 3a, b that the results obtained from the proposed first-order peridynamic homogenization theory correlate well with the results from [64,69].
To quantitatively measure the degree of anisotropy that may be induced on the elastic property of material due to crack growth (22) is utilized to compute the elastic anisotropy index A of the effective material tensor. Figure 4 shows the evolution of the elastic anisotropic indices for the case of a single crack and two coalescing cracks. As can be seen from the figure, the propagation of cracks results in increasing anisotropy. Comparison of anisotropic indices from the two cases shows an almost identical growth trend except for the jump corresponding to crack length close to when the two coalescing cracks were about to merge. This jump in elastic anisotropy index correlates well with the jump captured in the evolution of the C 22 component of the effective material stiffness tensor as shown in Fig. 3b.

Damage evolution due to randomly distributed microcracks
In this second example, we study the evolution of damage due to randomly distributed microcracks of equal length L . Two case studies would be investigated. In the first case, the RVE consists of microcracks that are aligned horizontally (Fig. 5a) while in the second case, the microcracks in the RVE are randomly aligned (Fig. 5b). In both case studies, the material of the RVE is taken to consist of an isotropic matrix having an elastic modulus E = 2 × 10 9 Pa and Poisson's ratio v = 0.3.
To obtain the effective properties of the material with evolving microstructure as characterized by crack growth, the implementation strategy highlighted in Sect. 3 was executed. The evolution of damage due to randomly distributed cracks is presented in Fig. 6a for cracks that are horizontally aligned and in Fig. 6b for cracks that are randomly oriented. Each curve represents the evolution of a normalized damage parameter n which measures the softening of the elastic modulus and is given as where E * n is the effective elastic modulus in the coordinate direction n = 1, 2 , such that n = 1 corresponds with the horizontal coordinate axis and n = 2 corresponds to the vertical coordinate axis. E is the elastic modulus of the undamaged matrix material.
In the case of horizontally aligned cracks, the damage parameter 1 remained largely unchanged as the damage process progresses as shown in Fig. 6a, which suggests that the material essentially retained its mechanical properties in the n = 1 direction. However, in the n = 2 coordinate direction, the evolution of the damage parameter 2 shows progressive softening of the material as the damage process progresses.
In the second scenario where the randomly distributed microcracks are given random orientations, the results of the computational homogenization as presented in Fig. 6b indicate that while the damage due to increasing prevalence of microcracks caused the material to soften, the effective material behaviour did not, however, deviate much from its original isotropic state.
It is inferred from Fig. 6a that the prevalence of horizontally oriented, randomly distributed microcracks transformed a hitherto isotropic material into an effective anisotropic material. To quantitatively measure the degree of anisotropy induced, (22) is employed to evaluate the elastic anisotropy index A . Figure 7 shows the elastic anisotropy indices for the two cases of horizontally aligned and randomly aligned cracks as the cracks accumulate. Figure 7 indicates that while the accumulation of horizontal cracks has the effect of increasing the anisotropy of the material, the accumulation of randomly aligned cracks essentially did not induce much anisotropy in comparison.
For the purpose of model validation, the evolution of effective damage obtained using the PDCHT was compared with previous results from OSBPDHT [64]. The results shown in Fig. 6a, b show agreement in the overall trend of the results from both methods. However, there are noticeable differences in the details of the results because although the same crack densities were used in both studies, the randomness in the distribution and orientation of microcracks in both studies made it impossible to exactly reproduce the distribution of the microcracks and their orientations as in [64]. Microcracks have been shown to interact with one another when their distance apart is less than some minimum value as demonstrated in the results presented in Fig. 3b. Since each instance of randomly generated microcracks gives rise to different arrangement and hence different interaction between microcracks, it is expected that analysis of two different distribution will yield different results. It has also been demonstrated in the results presented in Fig. 6 that the orientation of the microcracks has significant influence on the effective response. Since two instances of randomly generated cracks with arbitrary orientation give rise essentially to two different microstructures, it follows that the variation in the result of analysis from these instances of the microstructure is an expected consequence.

Conclusion
In this contribution, the recently proposed PDCHT was utilized to study the evolution of effective material response of damaged media. The PDCHT is an energy-based homogenization procedure which ensures constitutive equivalence between two structurally dissimilar materials by enforcing a strain energy equivalence between them. This computational homogenization was performed in a quasi-static analysis framework in which different loading conditions are assumed to yield progressive damage either in the form of propagating crack or increasing incidence of microcracks.
Two numerical examples were solved. The first example allows us to study the effect of propagating crack and crack interaction on the effective material coefficients. It was demonstrated that the presence of crack changes material behaviour in this case from isotropic to mild and strongly anisotropic. As cracks propagate, the material exhibits softening in one or both directions.
In the second example, we studied the damage evolution of a material due to proliferation of microcracks. Two case studies were investigated. The first case is when the Fig. 6 Evolution of damage due to random cracks. a Horizontally aligned, b randomly aligned Fig. 7 Elastic anisotropy index for cracks aligned horizontally and randomly microcracks are randomly distributed but horizontally aligned while the second case is where the randomly distributed microcracks are given arbitrary orientations. Effective material coefficients obtained from the first scenario show a strongly anisotropic effective behaviour while in the second case, the effective behaviour essentially remained isotropic. However, in both cases, the material exhibited softening with increasing damage accumulation.
A further development of these studies is expected. The foundation laid here can be utilized to study phenomena that occur out of interaction of microcracks with macro-cracks. Example of such phenomena is crack shielding and amplification of macrocracks due to interaction with microcracks. In this framework, the damage due to microcracks can be homogenized so that their interaction with the macrocrack would be applied in an average sense.