Effective Properties of Composites with Periodic Random Packing of Ellipsoids

The aim of this paper is to evaluate the effective properties of composite materials with periodic random packing of ellipsoids of different volume fractions and aspect ratios. Therefore, we employ computational homogenization. A very efficient MD-based method is applied to generate the periodic random packing of the ellipsoids. The method is applicable even for extremely high volume fractions up to 60%. The influences of the volume fraction and aspect ratio on the effective properties of the composite materials are studied in several numerical examples.


Introduction
Composite materials are widely used in engineering applications. Extraction of the mechanical properties by experiments is often expensive, time consuming and sometimes unfeasible. Therefore, it is important to develop modelling approaches, such as computational homogenization, to extract the mechanical properties of composite materials. The mechanical properties of composite materials are strongly determined by their microstructure, e.g., the shape, size and distribution of the inclusions in the matrix. Generating complex microstructures for RVEs (representative volume elements), which are commonly used in computational homogenization, remains a challenge.
Computational multiscale methods, which are commonly used for extracting the mechanical properties of composites, can be classified into hierarchical, semi-concurrent and concurrent methods [1]. A classical hierarchical multiscale method transfers information from the fine scale to the coarse scale only. It is particularly efficient for linear materials and becomes expensive for non-linear material behavior, since all possible deformation states have to be accounted for in order to extract the material response. The Cauchy-Born rule or the FE 2 approach [2,3] is a classical semi-concurrent method, which is also based on RVEs. However, in the FE 2 approach, information is passed back also from the coarse scale to the fine scale in the course of a simulation. Though semi-concurrent multiscale methods are computationally expensive, they account only for deformation states that actually occur in a simulation, and the constitutive model is determined on the fly. Concurrent multiscale methods also transfer information back from the coarse scale to the fine scale. However, concurrent multiscale methods do not exploit RVEs. Instead, the fine scale is directly embedded into the coarse scale and coupling via methods such as Lagrange multipliers, etc. There are numerous concurrent multiscale methods, such as the handshake domain method [4][5][6] or interface coupling methods [7].
where l micro , l macro are the length of the micro-and macro-scales, respectively. The second assumption is the periodicity of the microstructure. Two different scales x and y, associated respectively with the macroscale and microscale domains, are assumed. The relation between the microscale and macroscale spatial coordinates is: where x i (i = 1, 2, 3) denotes the component of macroscale coordinates; y i corresponds to the microscale coordinates; ε 1 is a scale factor between the macroscale and the microscale.
Consider a macroscale domain Ω ε with boundary Γ ε , which has a microstructure associated with a microscale domain denoted as Θ. The governing equation of the linear elastic problem may be expressed as: The displacement and traction boundary conditions applied to the domain are: where the subscript i, j, m, n ∈ {1, 2, 3}, σ ε ij , e ε mn , u ε i , L ε ijmn are the components of stress, strain, displacement and constitutive tensor, respectively, b i is the body force vector, u i , t i are the prescribed displacement and traction, Γ ε u , Γ ε t are the displacement and traction boundary and n j is the outward normal vector of the boundary. The subscripts x i , y i refer to the spatial coordinates of the macroscale and microscale, respectively.

The Asymptotic Expansion of Displacement
Assuming the material responses at the macroscale model denoted as x and the microscale model denoted as y are related, the displacement, strain and stress fields can be approximated by using the asymptotic expansions [23] in terms of ε, the scale factor between the micro-and macro-scales, where u (r) i is the y-periodic function of the order r of the displacement field. Substituting the asymptotic expansion of the displacements in Equation (8) into the strain-displacement relations in Equation (5), regrouping the terms according to the Equation (9), then the strain terms can be derived as follows, By substituting the strain in Equation (9) into the constitutive Equation (4), we have the asymptotic expansion of the stress as: The asymptotic expansion of the equilibrium equation can be obtained by substituting the stress in Equation (10) into the equilibrium Equation (3) and regrouping the terms according to the order of ε: After obtaining the equilibrium equations and boundary conditions related to different orders of ε, the differential equations of the two-scale problems can be derived by solving Equations (15)- (17). After solving the equations at different scales, the multiscale problem is decomposed into: The macroscale problem: The microscale problem: 19) and the bridging between the two scales is given by: The computational homogenization method may be divided into four steps [42] as illustrated in Figure 1, which has been previously implemented with convergence studies in ABAQUS [23]: Step 1. Solve the microscale problem in Equation (19) with six unit right-hand side strain vectors, and compute the stress σ kl ij with the associated strain.
Step 3. Solve the macroscale problem in Equation (18) with the homogenized elastic tensor.
Step 4. Post-process the local stress of the critical RVE. Step 2. Compute the homogenized elastic tensor ijkl L by Equation (20).
Step 3. Solve the macroscale problem in Equation (18) with the homogenized elastic tensor.
Step 4. Post-process the local stress of the critical RVE.

RVE Generation Algorithm
The shapes of the inclusions in the microstructure are various, and they affect the effective properties of random two-phase materials [43]. It is hard to describe the shape precisely and improve the computational efficiency at the same time. The ellipsoids with different aspect ratios can describe a wide range of shapes; see Figure 2. The sphere and normal ellipsoid can simulate 'regular' reinforcements. Fiber reinforcements may be simulated by prolate ellipsoids. Joints in the rock mass, microcrack or two-dimensional reinforcement, such as graphene, are usually described as discs, so they can be simulated by the oblate ellipsoids. In this paper, the generation algorithm of periodic random packing of ellipsoids will be presented, which is a molecular dynamics (MD)-based method. All of the centroids of the N ellipsoids are randomly created in a cell of side L at the beginning. Each ellipsoid has a random orientation and initially a null volume. The ellipsoidal volumes are gradually increasing at a constant growth rate. The growth rates of the semi-principal axes (a0, b0, c0) are chosen so that a0/a = b0/b = c0/c, where (a, b, c) are the semi-principal axes' lengths of the ellipsoids. At the same time, each ellipsoid has a random linear and angular velocity. The ellipsoids are then set in the translational and rotational motion, and the volumes gradually increase. At each increment, two types of collisions need to be checked, which are the collision between two ellipsoids and the collision between an ellipsoid and a cell face. If the collision between two ellipsoids occurs, the linear and angular velocities of the involved ellipsoids are updated according to the law of the conservation

RVE Generation Algorithm
The shapes of the inclusions in the microstructure are various, and they affect the effective properties of random two-phase materials [43]. It is hard to describe the shape precisely and improve the computational efficiency at the same time. The ellipsoids with different aspect ratios can describe a wide range of shapes; see Figure 2. The sphere and normal ellipsoid can simulate 'regular' reinforcements. Fiber reinforcements may be simulated by prolate ellipsoids. Joints in the rock mass, microcrack or two-dimensional reinforcement, such as graphene, are usually described as discs, so they can be simulated by the oblate ellipsoids. In this paper, the generation algorithm of periodic random packing of ellipsoids will be presented, which is a molecular dynamics (MD)-based method. Step 2. Compute the homogenized elastic tensor ijkl L by Equation (20).
Step 3. Solve the macroscale problem in Equation (18) with the homogenized elastic tensor.
Step 4. Post-process the local stress of the critical RVE.

RVE Generation Algorithm
The shapes of the inclusions in the microstructure are various, and they affect the effective properties of random two-phase materials [43]. It is hard to describe the shape precisely and improve the computational efficiency at the same time. The ellipsoids with different aspect ratios can describe a wide range of shapes; see Figure 2. The sphere and normal ellipsoid can simulate 'regular' reinforcements. Fiber reinforcements may be simulated by prolate ellipsoids. Joints in the rock mass, microcrack or two-dimensional reinforcement, such as graphene, are usually described as discs, so they can be simulated by the oblate ellipsoids. In this paper, the generation algorithm of periodic random packing of ellipsoids will be presented, which is a molecular dynamics (MD)-based method. All of the centroids of the N ellipsoids are randomly created in a cell of side L at the beginning. Each ellipsoid has a random orientation and initially a null volume. The ellipsoidal volumes are gradually increasing at a constant growth rate. The growth rates of the semi-principal axes (a0, b0, c0) are chosen so that a0/a = b0/b = c0/c, where (a, b, c) are the semi-principal axes' lengths of the ellipsoids. At the same time, each ellipsoid has a random linear and angular velocity. The ellipsoids are then set in the translational and rotational motion, and the volumes gradually increase. At each increment, two types of collisions need to be checked, which are the collision between two ellipsoids and the collision between an ellipsoid and a cell face. If the collision between two ellipsoids occurs, the linear and angular velocities of the involved ellipsoids are updated according to the law of the conservation All of the centroids of the N ellipsoids are randomly created in a cell of side L at the beginning. Each ellipsoid has a random orientation and initially a null volume. The ellipsoidal volumes are gradually increasing at a constant growth rate. The growth rates of the semi-principal axes (a 0 , b 0 , c 0 ) are chosen so that a 0 /a = b 0 /b = c 0 /c, where (a, b, c) are the semi-principal axes' lengths of the ellipsoids. At the same time, each ellipsoid has a random linear and angular velocity. The ellipsoids are then set in the translational and rotational motion, and the volumes gradually increase. At each increment, two types of collisions need to be checked, which are the collision between two ellipsoids and the collision between an ellipsoid and a cell face. If the collision between two ellipsoids occurs, the linear and angular velocities of the involved ellipsoids are updated according to the law of the conservation of momentum. If the collision between an ellipsoid and a cell face occurs, the periodic images of the involved ellipsoid are created on the opposite face, line or vertex. The process stops until the volume fraction V f is reached. The flowchart of the MD-based generation procedure of the periodic random packing of ellipsoids is shown in Figure 3.
Materials 2017, 10, 112 6 of 18 of momentum. If the collision between an ellipsoid and a cell face occurs, the periodic images of the involved ellipsoid are created on the opposite face, line or vertex. The process stops until the volume fraction Vf is reached. The flowchart of the MD-based generation procedure of the periodic random packing of ellipsoids is shown in Figure 3.

Representation of an Ellipsoid
An ellipsoid can be represented by the centroid vector r and the orientation denoted by two Euler angles, the dip angle θ and dip direction φ. The quaternion [44] is more convenient and stable

Representation of an Ellipsoid
An ellipsoid can be represented by the centroid vector r and the orientation denoted by two Euler angles, the dip angle θ and dip direction ϕ. The quaternion [44] is more convenient and stable than the two Euler angles in numerical calculations. The quaternion of the ellipsoid consists of a scalar and a vector. The vector is the rotation axis of the ellipsoid, and the scalar is the rotation angle. The quaternion of the ellipsoid i at time t can be written as: where: At the end of the procedure, if ellipsoid i at time t has the quaternion as Equation (21), the two Euler angles θ and ϕ are given by: where ψ i kt (k = 1, 2, 3) is the k-th term of the vector ψ i t . Assume that the global coordinate system is Ox 1 x 2 x 3 , and the local coordinate system O'x' 1 x' 2 x' 3 is aligned along the principal axes of ellipsoid i. The vector r i t denotes the position of the centroid O' at time t. The equation of ellipsoid i at time t in its local coordinate system can be written as: where: is a diagonal matrix where a t i , b t i , c t i denote the length of the semi-principal axes of ellipsoid i at time t; x = (x 1 , x 2 , x 3 , ω ) T are the homogeneous coordinates of a point, which are more convenient in the computation than the Cartesian coordinates. The equivalent Cartesian coordinates of x' are In homogeneous coordinates, the transition from the local coordinates system to the global coordinates system is directly combined by the rotation and translation, i.e., where the transition matrix M i t can be obtained from the quaternion of ellipsoid i at time t. R i t is the rotation matrix given by: where I is the identity matrix and S i t is an antisymmetric matrix given by [44]: Substituting Equations (28)- (30) into Equation (26), the equation of ellipsoid i at time t in the global system becomes: where:

Moving and Growing an Ellipsoid
The configuration of ellipsoid i at time t is determined by the position vector r i t , the quaternion q i t and the semi-principal axes' lengths (a i t , b i t , c i t ). The ellipsoid i has a linear velocity ν i and rotation velocity ω i . The growth rates of the ellipsoid i are (a i 0 , b i 0 , c i 0 ). The configuration of ellipsoid i at time t + ∆t can be derived according to the moving and growing rates of the ellipsoid i. The semi-principal axes' lengths at time t + ∆t are: The new position of the centroid at time t + ∆t can be derived by the linear velocity, given by: From time t to time t + ∆t, the ellipsoid i has been rotated by an angle ω i ∆t about the unit vector ω i ω i . This motion can be expressed in terms of a quaternion q i ∆t , given by: The quaternion at time t + ∆t can be obtained by combining the quaternion q i t and q i ∆t , such that [44]: where: By knowing the position vector r i t+∆t , the quaternion q i t+∆t and the semi-principal axes' lengths (a i t+∆t , b i t+∆t , c i t+∆t ), the coefficient matrix A i t+∆t can be deduced from Equation (27), and the transition matrix M i t+∆t is obtained from Equations (29) and (30). Finally, the equation of ellipsoid i at time t + ∆t is given by Equations (31) and (32).

Collision Time between Two Ellipsoids
The algebraic condition for the configuration relationship of two ellipsoids i and j can be represented by the roots of the characteristic equation of the two ellipsoids [45]. The characteristic equation at time t can be written as: where A i t and A j t are the coefficient matrix of ellipsoids i and j, respectively. Equation (39) is a fourth order polynomial of λ and has four roots. The relationship between the four roots and the configuration of the two ellipsoids is shown in Table 1 Table 1 is satisfied. This is an optimization problem, which can be solved with the bisection method and the secant method. Once the collision time t c is obtained, the contact point x c of the two ellipsoid can be solved from the following equation [46]: where λ 0 is the double positive root of the characteristic Equation (39) at time t + t c .

Collision Time between an Ellipsoid and Cell Faces
The point coordinates x in the equation x T A i t x = 0 of the ellipsoid i are a quaternion, but the point coordinates in the equation of a cell face are a vector. Therefore, the homogeneous coordinates need to be transformed to Cartesian coordinates. The coefficient matrix A i t can be rewritten as: where B i t is a 3 × 3 matrix; d i t is a 3 × 1 vector; F i t is a scalar. The equation of ellipsoid i becomes: The intersection of ellipsoid i with a plane x k = β (k ∈ {1, 2, 3}; β ∈ {0, L}) is an ellipse, a point or the empty set, which is shown in Figure 4. The intersection type depends on the spatial relationship of the plane with Points A and B, where the normal of the ellipsoid is parallel with the normal of the plane. The normal of the ellipsoid can be derived from the gradient of the ellipsoid equation, such that: where n(x) is the normal of the ellipsoid at point x. The intersection of ellipsoid i with a plane xk = β (k ∈ {1, 2, 3}; β ∈ {0, L}) is an ellipse, a point or the empty set, which is shown in Figure 4. The intersection type depends on the spatial relationship of the plane with Points A and B, where the normal of the ellipsoid is parallel with the normal of the plane. The normal of the ellipsoid can be derived from the gradient of the ellipsoid equation, such that: (43) where ( ) n x is the normal of the ellipsoid at point x.
Since the normals nA and nB are parallel to the normal of the plane nf, the points xA and xB can be solved from the following equation: =0 n x n x+d n (44) Assuming that xAk < xBk, the intersection relationship of the ellipsoid i with the plane xk = β is shown in Table 2. The collision time between the ellipsoids i with a plane is the time t + ts when Condition 3 in Table 2 is satisfied. This is an optimization problem that can be solved with the bisection method and the secant method.

Cases
Conditions Ellipsoids Configuration Since the normals n A and n B are parallel to the normal of the plane n f , the points x A and x B can be solved from the following equation: Assuming that x Ak < x Bk , the intersection relationship of the ellipsoid i with the plane x k = β is shown in Table 2. The collision time between the ellipsoids i with a plane is the time t + t s when Condition 3 in Table 2 is satisfied. This is an optimization problem that can be solved with the bisection method and the secant method. Table 2. Conditions for the intersection relationship of an ellipsoid with a plane.

Cases
Conditions Ellipsoids Configuration

Update the Velocities
Assuming that the linear and angular velocities of the ellipsoids i and j before the collision are (ν i− , ω i− ) and (ν j− , ω j− ), the linear (ν i+ , ν j+ ) and angular (ω i+ , ω j+ ) velocities after the collision can be solved according to the law of the conservation of momentum. The friction and the energy loss are assumed negligible. An orthonormal basis (n c , t 1 , t 2 ) is defined at the contact point x c , which is shown in Figure 5. n c is the normal of ellipsoid i at contact point x c , which is defined as Equation (43). The unit vector along n c is given by: Materials 2017, 10, 112 11 of 18 x + d n x + d (45) The unit vectors t1 and t2 are defined such that (nc, t1, t2) form a orthonormal basis, given by: (47)  The linear momentum along the vectors (nc, t1, t2) and the angular momentum about the contact point xc are conserved during the collision. The conservation of the linear momentum allows us to write the following equations: The unit vectors t 1 and t 2 are defined such that (n c , t 1 , t 2 ) form a orthonormal basis, given by: The linear momentum along the vectors (n c , t 1 , t 2 ) and the angular momentum about the contact point x c are conserved during the collision. The conservation of the linear momentum allows us to write the following equations: where r ∈ {i, j}; m r t is the mass of the ellipsoid r at time t with a unit density, which is given by: Considering the effect of the ellipsoids growing rate, the closing velocity of the two ellipsoids along vector n c must satisfy the following equation in order to avoid the ellipsoids from separating after the collision.
where ν c− and ν c+ are the closing velocities from ellipsoid i to j before and after the collision, respectively, given by: When an ellipsoid intersects with the cell faces at time ts, the periodic images depend on the number of faces that intersect the ellipsoid. When the number of faces is 1, 2, 3, the number of the periodic images is 1, 4, 7, respectively, which is shown in Figure 6. The periodic images have the same quaternions, velocities and growing rates as the ellipsoid intersecting with the cell faces. The centroid of the periodic images has an offset {−L, L}, which depends on the face for which the periodic images appear.

RVE Samples
Four RVE samples of different aspect ratios and volume fractions are generated with the MD-based method presented above, which is shown in Figure 7.

RVE Samples
Four RVE samples of different aspect ratios and volume fractions are generated with the MD-based method presented above, which is shown in Figure 7. Figure 7a-d shows the random packing of the sphere, normal ellipsoid, prolate ellipsoid and oblate ellipsoid, respectively. The aspect ratios are denoted as R 1 = a/b and R 1 = a/c, where semi-principal axis a is the axis of revolution. The number of inclusions denoted as N is the equivalent of the size of inclusions, when the size of the cell is fixed. During the numerical experiments, we tested increasingly larger samples of the material, while keeping the size of the cell fixed. It was found that the effective properties stabilized when approximately 20 inclusions were used in a sample, which will be discussed in Section 4; see Figure 8. Therefore, the samples in Figure 7 satisfy the assumptions of computational homogenization. is fixed. During the numerical experiments, we tested increasingly larger samples of the material, while keeping the size of the cell fixed. It was found that the effective properties stabilized when approximately 20 inclusions were used in a sample, which will be discussed in Section 4; see Figure 8. Therefore, the samples in Figure 7 satisfy the assumptions of computational homogenization.

Effective Properties of Composite Materials
We now want to determine the effective elastic properties of the composite material consisting of the periodic random packing of ellipsoidal inclusions embedded in an isotropic elastic matrix. The size, aspect ratio and volume fraction of the ellipsoids will affect the effective properties of the composite material. These parameters are controlled with the RVE generation algorithm. The effective properties of the composite are calculated by computational homogenization as presented in Section 2. The results depend on the mechanical parameters of the inclusions and matrix. The Lame constants of the inclusions and matrix are (λ 1 , µ 1 ) = (1, 1) and (λ 2 , µ 2 ) = (10, 10), respectively. All of the RVEs in this section are assumed as the unit cell and having the normalized Lame constants; thus, the output is the normalized value of the homogenized parameters. The normalized bulk modulus K K 0 and shear modulus G G 0 are chosen as the output, where K, G are the homogenized modulus and K 0 , G 0 are the modulus of the matrix.

Influence of the RVE Size
Firstly, the size of the representative volume element (RVE) needs to be determined. The size is acceptable when the effective properties of composite materials stabilize if the size is increasing. Since the cells are the unit in all of the RVEs, the RVE size is equivalent to determine the minimal number of ellipsoids, which is sufficiently large to include a sampling of the microstructure and is homogeneous in the sense of macroscopic effective properties. Two cases are considered to determine the minimal number of the ellipsoids. In the first case, the aspect ratios are R 1 = R 2 = 2, and the volume fraction is V f = 30%. In the second case, the aspect ratios are R 1 = R 2 = 10, and the volume fraction is V f = 10%. Five numbers N (5, 10, 20, 40, 60) are studied for each case. The normalized bulk and shear modulus are shown in Figure 8, which shows that the effective properties stabilize when the number of inclusions exceeds 20. The number of 20 is valid when the aspect ratio R 1 , R 2 ≤ 10. Indeed, it is necessary to generate more inclusions to obtain a uniform distribution of ellipsoid orientations when the aspect ratios are higher. This is shown in Figure 8b. The inclusion number of 20 will be considered in all of the following computations.

Influence of the Aspect Ratio
After the minimal number of the inclusions is determined, the influence of the aspect ratio on the effective properties of composite materials will be considered. As we know, the reinforcement effect mainly depends on the volume fraction when the macroscopic properties of the composite materials are isotropic. However, the aspect ratio will affect the reinforcement effect with different volume fractions. Two series of tests will be considered to study the influence of the aspect ratio on the reinforcement effect. Two cases are considered in the first series of the tests. In the first case, the aspect ratios are R 1 = R 2 = 2, and the inclusion number is N = 20. In the second case, the aspect ratios are R 1 = R 2 = 5, and the inclusion number is N = 20. Four volume fractions V f (5%, 10%, 15%, 20%, 25%, 30%) are studied for each case. The normalized bulk and shear modulus of the first series of test are shown in Figure 9. Three cases are also considered in the second series of tests. In the first case, the volume fraction is V f = 10%. In the second case, the volume fraction is V f = 20%. In the third case, the volume fraction is V f = 30%. Since the max volume fraction for the aspect ratios R 1 = R 2 > 10 cannot reach 30%, so the aspect ratios R 1 = R 2 (2,3,4,5,6) are studied for the volume fractions 10% and 30%. Additionally, the aspect ratios R 1 = R 2 (8,10,12,14,16) are studied for the volume fractions 10% and 20%. For the aspect ratios R 1 = R 2 (2,3,4,5,6) and (8,10,12,14,16), the inclusion numbers are N = 20 and N = 60, respectively. The normalized bulk and shear modulus of the second series of tests are shown in Figure 10.  Figure 9b. However, the composite materials are better reinforced with a higher aspect ratio with an increasing aspect ratio from 1-5 particularly for Vf = 30%. The influence of the aspect ratio on elastic properties has been studied in GP Tandon's paper [31]. It is observed that the Young's modulus and shear modulus increase with the increasing aspect ratio, whereas the bulk modulus decreases with it. These phenomena occur in the plane strain problem. The influence of the aspect ratio on elastic properties in the 3D microstructure is different. From Figure 10, it can be clearly observed that with the increase of the aspect ratio of the inclusions, there are no significant variations in the normalized shear modulus. On the other hand, the normalized bulk modulus increases slightly with increasing aspect ratio from 1-10. After 10, the influence of the aspect ratio gradually decreases. It also can be observed that when the volume fraction decreases from Vf = 30% to Vf = 20%, the normalized bulk and shear modulus significantly decreased. Generally speaking, the influence of the volume fraction on the elastic properties is greater than that of the aspect ratios.
(a) (b) Figure 9. The variation of normalized modulus with the volume fraction for the aspect ratio being fixed: (a) bulk modulus; (b) shear modulus.
(a) (b) Figure 10. The variation of normalized modulus with the aspect ratio for the volume fraction being fixed: (a) bulk modulus; (b) shear modulus.

Conclusions
In this paper, a computational homogenization model has been developed for the evaluation of the effective properties of the composite materials. The general steps of the computational homogenization model have been implemented in ABAQUS. The composite materials have the microstructure of the periodic random packing of ellipsoid inclusions, which are generated with an  Figure 9b. However, the composite materials are better reinforced with a higher aspect ratio with an increasing aspect ratio from 1-5 particularly for Vf = 30%. The influence of the aspect ratio on elastic properties has been studied in GP Tandon's paper [31]. It is observed that the Young's modulus and shear modulus increase with the increasing aspect ratio, whereas the bulk modulus decreases with it. These phenomena occur in the plane strain problem. The influence of the aspect ratio on elastic properties in the 3D microstructure is different. From Figure 10, it can be clearly observed that with the increase of the aspect ratio of the inclusions, there are no significant variations in the normalized shear modulus. On the other hand, the normalized bulk modulus increases slightly with increasing aspect ratio from 1-10. After 10, the influence of the aspect ratio gradually decreases. It also can be observed that when the volume fraction decreases from Vf = 30% to Vf = 20%, the normalized bulk and shear modulus significantly decreased. Generally speaking, the influence of the volume fraction on the elastic properties is greater than that of the aspect ratios.
(a) (b) Figure 9. The variation of normalized modulus with the volume fraction for the aspect ratio being fixed: (a) bulk modulus; (b) shear modulus.
(a) (b) Figure 10. The variation of normalized modulus with the aspect ratio for the volume fraction being fixed: (a) bulk modulus; (b) shear modulus.

Conclusions
In this paper, a computational homogenization model has been developed for the evaluation of the effective properties of the composite materials. The general steps of the computational homogenization model have been implemented in ABAQUS. The composite materials have the microstructure of the periodic random packing of ellipsoid inclusions, which are generated with an It is shown in Figure 9 that the dependence of the effective properties of the composite materials on the volume fraction of the inclusions for different aspect ratios is similar. The normalized bulk and shear modulus increase with the increasing volume fraction from 5%-30%. In Figure 9a, it is shown that the normalized bulk modulus is slightly higher under the higher aspect ratios R 1 = R 2 = 5. However, there are no significant variations in the normalized shear modulus with different aspect ratios in Figure 9b. However, the composite materials are better reinforced with a higher aspect ratio with an increasing aspect ratio from 1-5 particularly for V f = 30%. The influence of the aspect ratio on elastic properties has been studied in GP Tandon's paper [31]. It is observed that the Young's modulus and shear modulus increase with the increasing aspect ratio, whereas the bulk modulus decreases with it. These phenomena occur in the plane strain problem. The influence of the aspect ratio on elastic properties in the 3D microstructure is different. From Figure 10, it can be clearly observed that with the increase of the aspect ratio of the inclusions, there are no significant variations in the normalized shear modulus. On the other hand, the normalized bulk modulus increases slightly with increasing aspect ratio from 1-10. After 10, the influence of the aspect ratio gradually decreases. It also can be observed that when the volume fraction decreases from V f = 30% to V f = 20%, the normalized bulk and shear modulus significantly decreased. Generally speaking, the influence of the volume fraction on the elastic properties is greater than that of the aspect ratios.

Conclusions
In this paper, a computational homogenization model has been developed for the evaluation of the effective properties of the composite materials. The general steps of the computational homogenization model have been implemented in ABAQUS. The composite materials have the microstructure of the periodic random packing of ellipsoid inclusions, which are generated with an MD-based method. Different shapes of inclusions are generated in the RVE samples, which show that the MD-based method can generate various kinds of reinforcements and a high volume fraction of inclusions. The influence of the size, aspect ratio and volume fraction of the inclusions on mechanical properties has been studied. These studies showed that the effective properties of composite materials will be stable when the number of the inclusions is high enough to ensure the isotropy of the materials. There is no significant influence on the normalized shear modulus with the increase of the aspect ratio when the volume fraction is fixed. However, for the normalized bulk modulus, the reinforcement effect varies considerably for a lower aspect ratio and stabilizes for a higher aspect ratio. It shows that the MD-based method and the computational homogenization method are efficient tools to evaluate the effective properties of composite materials. In future work, it is desirable to apply the MD-based method to generating the microstructure of rock mass joints and microcracks. Not only modeling the elastic properties, the developed modeling scheme can be extended for the modeling of other effective macroscopic properties, such as effective thermal conductivity and permeability [47].