3-Dimensional Bond-Based Peridynamic Representative Volume Element Homogenization

In this study, a 3-dimensional (3D) implementation of representative volume element homogenization using bond-based peridynamic formulation is presented. Periodic boundary condition is established by coupling the displacements of periodic point pairs. Homogenized (effective) material properties are obtained based on peridynamic displacement gradient tensor. The current approach is validated by considering a composite material without defects and comparing homogenized properties with results obtained from another homogenization approach. Next, the capability of the current approach is demonstrated by considering randomly generated cracks with arbitrary orientation and location. It can be concluded that the current approach can be an alternative approach to obtain 3-dimensional homogenized material properties for heterogeneous materials with defects.


INTRODUCTION
To perform structural analysis at macroscopic scale, homogenized properties are utilized to approximately represent the material behavior. Although this approach usually gives reasonably accurate results, microstructural effects due to voids and microcracks can have significant influence on the macroscopic material behavior, especially fracture characteristics. Moreover, recently introduced microstructured materials as a result of advances in 3D printing and similar technologies require the consideration of the microstructural effects even further. To treat each individual microstructural effect as part of the analysis is usually intractable and homogenization approach is necessary which can also consider the effect of voids and microcracks.
There are various homogenization approaches available in the literature including rules of mixtures [1]. Homogenization process is usually done by performing volume averaging. Reuss and Voigt approaches provide the lower and upper limits of homogenized (effective) material properties [2]. In addition, a combination of finite element based representative volume element (RVE) has been developed. More-over, asymptotic homogenization theory was introduced which is based on a two-scale formulation [3]. Although these approaches are effective for homogenization purposes, they can encounter difficulties to deal with defects which cause discontinuities in the formulation. Spatial derivatives used as part of these homogenization approaches are not defined along crack surfaces. In such cases, finite element analysis can be a suitable option. However, special elements should be introduced to define the singularity at the crack tip [4]. Moreover, suitable crack propagation criteria should be introduced. Worth to mention that asymptotic homogenization, when used in combination with phase field approach, could potentially yield higher order parameters when properly applied. As an alternative, newly developed peridynamic formulation [5] can be more suitable since equations of motion of peridynamics are in the form of integrodifferential equations and do not contain any spatial derivatives. Definition of failure is also pretty straightforward with respect to the ones used in traditional approaches. Peridynamics is suitable to analyze both isotropic, anisotropic [6] and polycrystalline [7] materials. It is possible to perform both static and dynamic fracture analysis [8]. Any number of cracks can occur in the material and it is possible to investigate the interaction of microcracks and macrocracks [9,10]. Various peridynamic beam, plate and shell formulations are also available in the literature [11][12][13]. Homogenization can also be done within peridynamic framework [14]. In addition to structural analysis, peridynamics can also be used for the analysis of other physical fields [15][16][17].
As an extension of our earlier two-dimensional studies, in this study, peridynamic representative volume element homogenization approach is extended to three-dimensional configurations within bond-based peridynamic framework. Therefore, it will be possible to perform homogenization by considered defects such as microstructure.
To validate the current approach, peridynamic homogenization is performed by considering a composite material and effected properties are calculated. Peridynamic results are compared against another homogenization approach. Finally, to demonstrate the capability of the current approach, randomly oriented cracks are introduced in representative volume element to investigate their influence on effective material properties.

Bond-Based Peridynamic Formulation
Peridynamics was introduced by Silling [5] by reformulating the equations of continuum mechanics. For the original formulation, named as bond-based peridynamics, the equations of motion of the material point x can be expressed as [6] ( ) ( , ) ( , , )d ( , ).
In Eq. (1), H k defines the domain of influence which is called as horizon. Therefore, material point x can interact with all material points x′ inside its horizon with a size of δ. Note that peridynamic equations reduce to classical continuum mechanics equations as the horizon size converges to zero, u and u′ are displacements of material points x and x′, respectively, ρ, , u  f and b represent density, acceleration, peridynamic force between two material points and body force, respectively.
It is usually not feasible to obtain analytical solution of the peridynamic equation given in Eq. (1). Therefore, numerical methods can be utilized based on meshless approach and Eq. (1) can be written in the discretized form as [6] ( ) ( , ) In Eq. (2), N represents the total number of material points denoted as j inside the horizon of the material point k, V j is the volume of the material point j (Fig. 1), f kj is the peridynamic force between material points k and j which is defined as According to both based peridynamics, peridynamic force that the material point k is acting on material point j and the peridynamic force that the material point j is acting on material point k are equal in magnitude and opposite in direction. Therefore, both the conservations of linear and angular momentum are automatically satisfied. Moreover, the peridynamic force is only dependent on the motion of interacting material points.
The material constant in bond based peridynamics c is called bond constant and can be represented in terms of elastic modulus of the material E as 4 12 .
E c   (4) Note that in bond-based peridynamics, the value of Poisson's ratio is constrained to the value of ν = 1/4. The stretch of a peridynamic interaction (bond) s is defined as  where y k and y j are new positions of the material points x k and x j , in the deformed configuration.
In peridynamic formulation, since material points can interact with each other, interacting material points k and j can correspond to different materials. In such cases, the bond constant of the interaction between material points k and j can be obtained by splitting the interaction into multiple pieces (Fig. 2a) and can be calculated as where M represents the number of pieces and d m is the length of each piece of the interaction. For each piece, the bond constant c (k, j; m) can be calculated by using Eq. (4).
If an interaction is passing through a crack surface, then this bond will be broken (Fig. 2b). To represent this, the parameter 0 kj  is introduced which is defined as 0 0, if interaction passes through a crack surface, 1, otherwise.

Representative Volume Element Formulation
In this section, 3-dimensional bond-based peridynamic representative volume element homogenization will be presented. To accomplish this, a representative volume element (RVE) will be defined to approximate the microscopic material properties of heterogeneous medium (microscale) as material properties of a statistically homogeneous medium (macroscale).
The macroscale will be denoted as (x) and the microscale will be denoted as (y). Only linear elastic deformations will be considered in this study. For a heterogeneous material at the microscale, the constitutive relationships can be expressed as [2] , ,       C S (8) where S and C represent compliance and stiffness tensors, respectively, and ε and σ represent microscopic strains and stresses, respectively.
Under the assumption of microstructure being periodic, a unit sized RVE can be utilized for the micromechanical analysis which corresponds to a material point in macroscale analysis. After the homogenization process, constitutive relationships can be represented at macroscale as [2] * * , ,       C S (9) where S * and C * represent homogenized (effective) compliance and stiffness tensors, respectively, and ε and σ represent macroscopic strains and stresses, respectively.
According to the homogenization theory, the macroscopic displacement u can be expressed in terms of the microscopic displacement field u as [2] 1 d , where V is the volume of the RVE. Microscopic displacement field u calculated as , where u is the volume averaged displacement, u  is the displacement fluctuation function, and the average strain vector ε is defined as The macroscopic stress and strain can be expressed in terms of microscopic stresses and strains via volume averaging process as [2] 1 A constant strain tensor ij  is applied on a fictitious boundary region in peridynamic analysis and periodic boundary region relationship are expressed as shown in Fig. 3. In Fig. 3, the outer fictitious regions are called as "child". The material points inside "child" regions do not have their horizon. However, they can be considered as a material point within the horizon of nonfictitious material points. The regions next to "child" regions are called "parent" regions. Each material point inside a "child" region is coupled with a material point inside a "parent" region. On the other hand, each material point inside a "parent" region can be linked to multiple material points in "child" regions. The displacement of a material point inside a "child" region depends on both the periodic condition and the displacement of the associated material point in "parent" region. In contrast, the displacement of a material point in a "parent" region depends on material points inside of its horizon regardless of these material points being inside the fictitious or non-fictitious regions. In this study, the size of the fictitious region is considered as larger than twice of the horizon size so that all material points in the non-fictitious region has complete horizon and the periodic condition is satisfied.
The periodic relationship between material points inside "child" and "parent" regions can be expressed as where superscripts p and c denote paired material points inside "parent" and "child" regions respectively (Fig. 3), and ij  is the macroscopic strain.
By using Eq. (11) and periodic boundary condition given in Eq. (14), the microscopic strain field ε can be obtained as .   ε ε ε  (15) If we substitute Eq. (13) into the constitutive relations given in Eq. (9), the relationship between microscopic stress and macroscopic strain can be obtained as (16) in which the macroscopic strain  is applied to the cell as boundary condition and can be defined as The microscopic stress tensor σ can be obtained from the displacement-based homogenization result. As given in Breitenfeld et al. [18], the peridynamic stress vector s at material point k can be obtained from displacement vector U as  (18) in which D is the isotropic elastic moduli matrix, K is the 3 × 3 shape tensor stored in a 6 × 9 matrix [18], and the definition for N is given in [18], which is a 9 × 3n matrix with n being the number of material points within k's horizon. The displacement vector U is assembled as [18] [ . . .
where u, v and w are the displacement components in x, y and z direction, respectively. Peridynamic stress tensor σ can then be assembled from peridynamic stress vectors s for every material point. Finally, the effective stiffness tensor C * can be assembled column by column from the peridynamic stress tensor σ and the macroscopic strain ε as  where the angle bracket denotes the volume average over cell domain.

NUMERICAL RESULTS
The current approach is validated by considering a boron aluminum composite cell as shown in Fig. 4. The homogenized material properties from peridynamics are compared with the results from variational asymptotic method for unit cell homogenization (VAMUCH) [19]. The constituent properties for the boron aluminum composite are given in Table 1. In numerical analysis, a three-dimensional unit cell is considered with a discretization of 60 × 60 × 60. The spherical horizon size δ is specified as 3 times of discretization size. Homogenized (effective) material properties are obtained by using Eqs. (16)-(20) and tabulated in Table 2. According to these results, a good agreement is achieved between results from peridynamic homogenization approach and VAMUCH approach.
The second example considers a three-dimensional aluminum cube contains arbitrarily oriented cracks as shown in Fig. 5. The material properties for aluminum used in this example are identical to the aluminum properties listed in Table 1.  The cell geometry is considered to have unit dimension. The effective strength of such a cell is influenced by parameters such as crack orientation, crack density and crack shape etc. A pseudocrack density parameter is defined as where N is the total number of cracks, A is the surface area of an individual crack and V is the volume of the cell. All cracks share the same dimension as shown in Fig. 6, where w = 0.1 m and h = 0.05 m. However, the location and spatial orientation of each crack are randomly generated. Rule is set so that no generated crack penetrates outside cell domain.
Each crack has two unique properties: a location vector x(x, y, z) pointing towards the center of the crack, and a spatial rotation revolving around the crack center denoted as quaternion q(a, b, c, w). The purpose of representing crack rotation in a quaternion form is to accurately and conveniently obtain the Cartesian coordinates of four crack corner points. Although common Euler rotation can be applied, it tends to suffer from issues such as Gimbal lock and accuracy loss. Unit quaternions on the other hand does not have those limitations. The location vector can be simply generated using a random number ge- Fig. 6. Random three-dimensional crack location and orientation. nerator. However, in order to generate randomly distributed quaternions in a uniform way, a hypersphere in four-dimensional space called Hopf fibration is created and chosen as selection space for quaternions. Once randomly generated vector x and quaternion q are applied to the crack, the Cartesian coordinates of corner points c 0 , c 1 , c 2 and c 3 are calculated. Crack is numerically represented by simply breaking peridynamic bonds as shown in Fig. 2b. To determine which bond meets the breaking criteria require intersection check between bond and cracks. In this study, two separate ray casting tests are performed for each pair of bond/crack candidate to determine if the bond intersects with a given crack. The crack as shown in Fig. 6 is divided into two triangles, triangle c 0 -c 1 -c 3 and c 1 -c 2 -c 3 . One ray casting test is performed for each triangle.
Six crack densities are considered in this example: ρ = 0.1, 0.2, 0.3, 0.4, 0.5 and 0.6. Crack amount N is rounded into its nearest integer if necessary. Five instances are generated for each crack density level to obtain average results for comparison. List of crack  densities and their corresponding crack amount is shown in Table 3. A total of 30 numerical cases are performed. Examples are solved with 20 × 20 × 20 peridynamic discretization. Spherical horizon with δ equivalent to 3 times of discretization size is chosen.
Numerical results are shown in Figs. 7 and 8 where E x , E y and E z are reduced effective elastic moduli in x, y and z directions, respectively. The reduction in elastic modulus is resulted from the influence of cracks presented within the cell. The impact is measured by comparing the reduced modulus E x , E y and E z against their original undamaged counterpart 0 , x E 0 y E and 0 . z E The ratio 0 E E is denoted as an elastic softening parameter ψ E . Similarly, G xy , G yz and G xz are reduced effective shear modulus. Material softening parameter ψ G is defined as G/G 0 , where G 0 is undamaged shear modulus. When material softening happens due to the presence of cracks, parameter ψ E and ψ G will reduce accordingly.
For the randomly generated arbitrary oriented microcracks, the effective material properties in each direction decrease at the same rate while the crack density increases.

CONCLUSIONS
In this study, a 3-dimensional bond-based peridynamic representative volume element homogenization method for heterogeneous materials with defects is presented. The current approach is verified by considering a boron-aluminum composite cell and comparing homogenized material results with results obtained from another homogenization method. The homogenization of a unit cell with randomly generated arbitrarily orientated microcracks with different crack densities is also performed. It is observed that by increasing crack density, effective stiffness reduces at the same rate in both x, y and z directions. Therefore, it can be concluded that the current approach can be a suitable alternative for homogenization of heterogeneous materials with defects. Future study can be done to improve current method with ordinary statebased PD formulation, which will provide more homogenization options in terms of material properties such as Poisson's ratio.