A Lattice Model for Elastic Particulate Composites

In the present article, a version of the lattice or spring network method is proposed to model the mechanical response of elastic particulate composites with a high volume fraction of spherical particles and with a much weaker matrix compared to the stiffness of the particles. The main subject of the article is the determination of the axial stiffnesses of the springs of the cell. A comparison of the mechanical response of a three-dimensional particulate composite cube obtained using the finite element method and the proposed methodology showed that the efficiency of the proposed methodology increases with an increasing volume fraction of the particles.


Introduction
The lattice or spring network method is applied widely in various areas of mechanics: to solve various problems of continuum mechanics, micromechanics, molecular dynamics, fracture mechanics, multiscale modelling, soft materials, and so on [1][2][3][4][5][6]. Apart from that, this method can also be applied to model different materials: metals, concrete, asphalt, ceramics, various composites, particulate solids, granular matter, and biomaterials [7][8][9][10]. In addition, the lattice method can be applied to elasticity and viscoelasticity problems [1], or used as an alternative to the finite element method [2].
The geometry of a cell that approximates a continuum, the model of the spring of the cell, and the determination of the required stiffness parameters of the spring are three of the most important issues in the lattice model. When springs are used as the connecting elements of the nodes of the cells, only one parameter of the springs, i.e., the axial stiffness, has to be determined in terms of the material properties. It is relatively easier to determine the stiffness parameters of the spring for homogeneous materials than for composites. One approach to model composites using the lattice model is to approximate the constituent parts of the composites-for example, the matrix and the particles-by distinct springs for each constituent part of the composite [8,11,12]. In [12], it is assumed that the composite, consisting of bonded spherical particles, is approximated by a network of one-dimensional springs connecting the centres of the interacting spherical particles. The total stiffness of the spring is calculated as the total stiffness of the three distinct springs concatenated in a series. Two springs, the first and third, correspond to the interacting particles, while the intermediate spring corresponds to the bond between the particles. A similar approach was proposed in [13]. In this article, the convex spherical surfaces of the particles and the concave spherical surfaces of the bond element were taken into account in the evaluation of the axial stiffness of the spring that connects the centres of the two interacting particles. In this article, the connecting spring was also modelled as three distinct springs concatenated in series.
In the present paper, the ideas of the evaluation of the stiffness of the spring that connects the centres of two interacting adjacent particles, presented in article [13], are extended by evaluating the lower and upper bounds of the stiffness of the spring. Two different approaches were applied to evaluate the stiffnesses: the conditional connecting element was divided into infinitely small prisms, and infinitely thin rings. The former provides the lower bound and the latter provides the upper bound. The obtained closed-form solutions of both bounds of the stiffnesses were verified using the finite element method by modelling a 3D elastic particulate composite with different particle volume fractions, different Poisson's ratios, and different ratios of the modulus of elasticity of the particles to the matrix.
A particulate composite consisting of particles embedded in the matrix is approximated by one-dimensional springs that connect the centres of the particles (Figure 1). The springs are characterized by their length R p + L c , where R p is the particle diameter and L c is the distance between the particles, and the axial stiffness is K s (Figure 1b). The following assumptions are valid for the particulate composite and its approximating spring network: • the particles, matrix and the connecting springs obey the linear elastic law, • the particles interact with the matrix by their entire surface, • the particles of the composite do not rotate, • the interface member connecting two adjacent particles is composed of the matrix and is cylindrical, • the diameters of the ends of the interface member are the same as of connected particles, • the interface members interact only with the two adjacent particles and interact with the particles only by the entire surface of the hemispheres, • the pin-connected springs connect the centres of the adjacent particles, • the spring is composed of the interface member and two hemispheres, and only carries an axial force.
All quantities related to the particles are denoted by the subscript p and the quantities related to the interface members are denoted by the subscript b. The particle is characterised by the radius R p , elasticity modulus E p , and Poisson's ratio v p . The interface member is characterised by the cylinder radius R b and the elasticity constants E m and v b , respectively. Two limits K s,I and K s,I I of the axial stiffness can be obtained by considering two different divisions of the connecting element: as parallel and sequentially connected springs ( Figure 1).

Governing Equations for the Lower Bound of the Stiffness
Let us consider the half of the connecting element comprised of a half of the particle and a half of the interface member to obtain the lower bound K s,I of the stiffness of the connecting element ( Figure 2). The entire hemisphere is denoted hereafter by Ω p and its circular basis by S p . The entire half of the interface member is denoted hereafter by Ω b and its circular basis by S b .
Let the circular bases S p and S b of the interface member and the hemisphere (Figure 2a) be divided into n infinitesimal rectangles S p,ξ and S b,ξ so that S i = n ξ=1 S i,ξ and n ξ=1 S i,ξ = ∅, where i ∈ {p, b}. Then, the areas ∆A i,ξ of the rectangles S i,ξ , i ∈ {p, b}, are such that A i = ∑ n ξ=1 ∆A i,ξ , where A i , i ∈ {p, b} is the total area of the circular basis S i . Let the hemisphere Ω p and the interface member Ω b (Figure 1a) be divided into infinitesimal parallel connected prisms ∆Ω i,ξ whose bases are the rectangles S i,ξ so that Ω i = n ξ=1 Ω i,ξ , i ∈ {p, b}; see Figure 2a. Then, the stiffness of the prisms Ω p,ξ and Ω b,ξ connected sequentially can be expressed as follows (in Figure 2a only half of the prisms are depicted): where D p and D m are the elastic constants of materials of the particle and interface member, respectively. In case of the uniaxial stress state of the particle and the interface member, i.e., when σ xx = σ yy = 0 of the conditional connecting element, see Figure 2, the elastic constant D p = E p and D m = E m . When it is assumed that, for the conditional connecting element xx = yy = 0, then Other cases of the elastic constants are also possible. In Equation (1), l p (x p,ξ , y p,ξ ) and l b (x b,ξ , y b,ξ ) are the half lengths of the prisms of the hemispheres Ω p,ξ and Ω b,ξ at the points (x p,ξ , y p,ξ ) ∈ S p and (x b,ξ , y b,ξ ) ∈ S b , respectively. The total stiffness of all prisms Ω ξ = Ω p,ξ Ω b,ξ , ξ ∈ {1, ..., n}, or the stiffness of the connecting element K s = ∑ n ξ=1 ∆k ξ . By letting ∆A i,ξ → 0, we obtain a limit K s,I = lim ∑ ξ ∆k ξ , where ∆k ξ is given in Equation (1). The obtained limit can be rewritten as a double integral over where the integrand T(x, y) . ( In Equation (3), the half lengths of the prisms l p (x, y) and The integral of Equation (2) can be expressed as an iterated double integral in the rectangular coordinate system in the polar coordinates where r = x 2 + y 2 . After a change of variable z = (R 2 p − r 2 ) in Equation (7), we obtain Integration of Equation (8) yields the stiffness of the connecting element as D p = D m where a = (L c /2 + R p )D p , and b = D m − D p . It should be noted that, when D m = D p , then K s,I by Equation (9) is undefined due to the division by zero, since b = D m − D p = 0. In this case, K s,I can be obtained from Equation (7) by letting D m = D p : It is easy to notice that the obtained limit of the stiffness K s,I in Equation (10) corresponds to the stiffness of the homogeneous cylinder, i.e., K s,I = AD/l = πR 2 p D/(2R p + L c ), where A is the area of the cross-section of the cylinder, and D is the elastic constant of the cylinder.

Governing Equations for the Upper Bound of the Stiffness
In this subsection, the upper limit K s,I I of the stiffness of the connecting element is obtained by dividing the hemisphere Ω p and the interface member Ω b (Figure 1b) into sequentially connected cylinders of the infinitesimal height ∆h. The stiffness of the sequentially connecting cylinders can be expressed as follows: where ∆k i is the stiffness of cylinder i of infinitesimal height ∆h In Equations (12) and (13), A p (z) and A b (z) are the areas of the cross-sections of the particle and the interface member, respectively, dependent on coordinate z: A p (z) = π(R 2 p − z 2 ) and A b (z) = πz 2 . Then, Equation (11) can be rewritten as follows: The limit of the sum lim (14) can be written as an integral Integration of Equation (15) depends on D p and D m where the arctan and arctanh are the inverse tangent and inverse hyperbolic tangent, respectively. When D p = D m , then it is not possible to use Equations (16) and (17), since D m − D p = 0. In this case, K s,I I can be obtained from Equation (15) by letting D m = D p Finally, the stiffness K s,I I is It is easy to see that, when D m = D p , Equations (10) and (19) are equal, i.e., K s,I = K s,I I .

Mathematical Analysis of the Stiffnesses
In this subsection, hereafter, an analysis of the obtained Equations (9), (10), and (19) is presented The stiffnesses K s,I and K s,I I are nonlinear with respect to D p , D m , L c and R p , where the nonlinearity of a function f is defined as where c = L c /2 + R p , and lim where b is given in the explanations of the notations below Equation (9), and lim It should be noted that K s,I I tends to K s,I with increasing L c (see Figure 3); therefore, at the big values of L c , the stiffnesses K s,I ≈ K s,I I (Figure 3). In Figure 3, it is also demonstrated that K s,I and K s,I I approach each other when D p → D m or D m → D p .

Numerical Validation of the Proposed Methodology
Two numerical validations of the proposed methodology are presented hereafter in two subsections. In the first validation, the stiffnesses K s,I and K s,I I are compared with the stiffnesses of the 3D FE models of the connecting elements. In the second one, the mechanical behaviour of the particulate composite cube approximated by the spring model (SM) and modelled by 3D FE is compared. The FE analysis was performed by ANSYS 12.

Stiffness of the Connecting Element
The obtained Equations (9), (10) and (19) of the stiffnesses K s,I and K s,I I were verified by a 3D FE analysis of the connecting element shown in Figure 4. Two cases were considered. In the first case, E p = 40 GPa, ν p = 0.0, E m ∈ {40 × 10 9 , 30 × 10 9 , 20 × 10 9 , 10 × 10 9 , 5 × 10 9 , 2.5 × 10 9 , 5 × 10 8 , 5 × 10 7 , 1.0 × 10 6 } Pa, ν m = 0.498. In the second case, E m = 40 GPa, ν m = 0.0, E p ∈ {40 × 10 9 , 30 × 10 9 , 20 × 10 9 , 10 × 10 9 , 5 × 10 9 , 2.5 × 10 9 , 5 × 10 8 , 5 × 10 7 , 1.0 × 10 6 } Pa, ν p = 0.498. For both cases, R p = R b = 5 × 10 −3 m, and L c ∈ {5 × 10 −6 , 5 × 10 −5 , 5 × 10 −4 , 2.5 × 10 −3 } m. These parameters of the particles and interface member were chosen so that it would be possible to verify the obtained equations at different moduli of elasticity E i , Poisson's ratios ν i , i ∈ {p, m}, and at the different distances between the surfaces of the particles L c . The sample under investigation and its FE model, a quarter of the sample consisting of 1.087620 × 10 6 nodes and 0.778623 × 10 6 elements, with mesh are depicted in Figure 4a,b, respectively. The boundary conditions were as follows: U z = 0 was applied to plane B, and U x = U y = 0 was applied to the centre line M (Figure 4a, dotted line) of the FE model. The FE model was discretised by the tetrahedron elements "SOLID187" of 10 nodes of six degrees of freedom ( Figure 4). The average size of the finite elements of the particles of the sample was 0.01R p . The volume of the interface material was conditionally divided into two regions. The contact region between spheres was meshed by fine mesh whose average size of the finite elements was 0.015R p , while the remaining volume was meshed by a coarser mesh of an average size 0.002R p (Figure 4b).
The stiffness obtained by the FE method, denoted hereafter as K s,FEM , was obtained by applying a displacement ∆l on free plane A ( Figure 4a) and was calculated as K s,FEM = F/∆l, where F is the total reaction force of plane B at the displacement ∆l.
The calculation results are shown in Figure 5. As we can see from Figure 5, when E m ≤ E p , then, in the majority of the examined cases, K s,FEM is closer to K s,I than to K s,I I , i.e., |K s,I − K s,FEM | < |K s,I I − K s,FEM |. However, when L c = 2.5 mm and E m ∈ {20, 30, 40} GPa, then K s,FEM is closer to K s,I I than to K s,I . In addition, from Figure 5, we can see that in the majority of cases, except from the case when L c = 2.5 mm and E m ∈ {20, 30, 40} GPa, the stiffness K s,FEM ∈ [K s,I , K s,I I ]. In Figure 5, it is clearly depicted that K s,I , K s,I I → K s,FEM when E m → 0.

Mechanical Behaviour of a Particulate Composite Cube
The developed SM was validated by comparing the mechanical responses of a 3D particulate composite (see Figure 6) obtained by the 3D FE method and by SM. The stiffness of the springs of SM was calculated by the developed formulas given in Equations (9) and (10).
Overall, 126 samples, which differ in the modulus of elasticity of matrix E m and the volume fraction of the particles φ p , were calculated. The properties of the 3D FE model and the SM are the following (see Figure 6): the diameters of the particles d p ∈ {10, 12, 13} mm and the corresponding volume fractions of the particles φ p ∈ {28.46, 51.45, 61.35} %; the dimensions of the cube are (see Figure 6): height h c = 40.2 mm, width b c = 46.42 mm and depth d c = 32.82 mm. The moduli of elasticity of the 3D FE model and the springs of the SM are the following: E m ∈ {40 × 10 9 , 30 × 10 9 , 20 × 10 9 , 10 × 10 9 , 1 × 10 9 , 1 × 10 8 , 1 × 10 7 , 1 × 10 6 } Pa for the matrix, and E p = 40 × 10 9 Pa for the particles. Poisson's ratio of the particles and matrix for the 3D FE model are ν p , ν m ∈ {0.0, 0.33}, and the Poisson's ratios for the particle and matrix are the same, i.e., ν p = ν m . The stiffnesses of the springs were calculated by Equations (9) and (10). The elastic constants of the springs for SM were taken as for the uniaxial stress state, i.e., D i = E i , i ∈ {p, m}. Therefore, Poisson's ratio does not affect the stiffnesses K s,I of the springs. The length of the connecting elements L c + d p and the distance between the particles L p of the samples depend on the volume fraction φ p : L c = 3.4 mm for φ p = 28.46% (d p = 10 mm), L c = 1.4 mm for φ p = 51.45% (d p = 12 mm) and L c = 0.54 mm for φ p = 61.35% (d p = 13 mm). It is determined that the stiffness K s,FEM is closer to K s,I than to K s,I I as L c is small enough. Therefore, the results calculated only with the stiffness K s,I are presented hereafter. The 3D FE model and SM of the composite consist of tetrahedron lattices.
To validate the proposed methodology, the vertical U y and horizontal U x displacements in the directions y and x were imposed to the top planes of the corresponding samples, see Figure 7, and the tensile F y,3D and shear F x,3D forces of the 3D FE model were compared with the corresponding tensile F y,1D and shear F x,1D forces of the FE model of SM. The boundary conditions for the 3D FE model and SM were as follows: the displacements of the bottom plane of the 3D FE model and SM were restricted fully, i.e., U x = U y = U z = 0, see Figure 7.  Two loading cases were applied to the specimens to validate the proposed methodology, see Table 1. For the sake of illustration, the shear displacements U x of the 3D FE model in the direction x subject to loading case LC2 are shown in Figure 8. The dependences of the tensile and shear forces F y,1D , F y,3D , F x,1D , and F x,3D , of the loading cases LC1 and LC2 on the ratio E m /E p calculated by the 3D FE method and by SM, when d p = 10 mm at different Poisson's ratios ν p , ν m ∈ {0.0, 0.33} and different particles diameters d p ∈ {13, 12, 10} mm in double logarithmic scales, are shown in Figure 9. It should be noticed that the calculated tensile and shear forces F y,1D and F x,1D of SM do not depend on Poisson's ratios ν p , ν m ∈ {0.0, 0.33}, since the calculations were performed as D i = E i , i ∈ {p, m}, i.e., D i does not depend on Poisson's ratio.
As we can see from Figure 9, the agreement between the results of F x,1D and F x,3D as well as between F y,1D and F y,3D seems very good in double logarithmic scale at various ratios E p /E b and ν ∈ {0.0, 0.33}. However, the relative ratios of the forces (F y,3D − F y,1D )/F y,3D and (F x,3D − F x,1D )/F x,3D can reveal the agreement between the results better. Figure 9. Dependences of the tensile forces F y,3D and F y,1D for (a); and the shear forces F x,3D and F x,1D for (b) of the loading cases LC1 and LC2 on the ratio E p /E m at different Poisson's rations ν p ∈ {0.0, 0.33} when particles' diameters d p = 10 mm.
The dependences of the relative ratios (F y,3D − F y,1D )/F y,3D and (F x,3D − F x,1D )/F x,3D of the loading cases LC1 and LC2 on the ratio E p /E m at different Poisson's rations ν p ∈ {0.0, 0.33, 0.495} and particle diameters d p ∈ {13, 12, 10} in semi-logarithmic scales are shown in Figure 10. The relative ratios of loading case LC1 are shown in Figure 10a,c,e while those for the loading case LC1 are shown in Figure 10b,d,f. The ratios shown in Figure 10a,b correspond to the case when d p = 13 mm, whereas, in (c) and (d), to the case when d p = 13 mm, and in (e) and (f) to the case when d p = 10 mm. As we can see from Figure 10, there is not any unique tendency for the ratios |(F x,3D − F x,1D )/F x,3D | and |(F y,3D − F y,1D )/F y,3D | dependent on E p /E m except for the fact that the variation of the relative ratios is smaller when E p /E m > 10 2 Pa. From Figure 10, we can also see that for the tension loading case LC1, when d p ∈ {12, 13} mm, the relative difference |(F y,3D − F y,1D )/F y,3D | is the biggest as ν p = ν m = 0.495. However, a similar conclusion is not valid for the shear loading case LC2. Only when d p = 13 mm and E p /E m ≥ 40, the ratio |(F y,3D − F y,1D )/F y,3D | is the biggest for ν m = ν p = 0.495. The value of the Poisson ratio 0.495 is an extreme case. It real life, for common materials, the Poisson ratio can be assumed as ν ∈ [0.15, 0.40]. The analysis showed that, for the calculated cases, the following limits are valid as d p ∈ {10, 12, 13} mm: for the loading case LC1 (F y,3D − F y,1D )/F y,3D ∈ [−0.16, 0.39] as ν = 0.0 and (F y,3D − F y,1D )/F y,3D ∈ [−0.11, 0.42] as ν = 0.33; while for the loading case LC2 (F x,3D − F x,1D )/F x,3D ∈ [−0.08, 0.27] as ν = 0.0, and (F x,3D − F x,1D )/F x,3D ∈ [−0.18, 0.12] as ν = 0.33. Since the width of the intervals of the relative ratios of the loading case LC2 are less than for the loading case LC1, then the proposed methodology is more accurate for LC2 than for LC1.
The figure clearly shows that, when D P = 13 mm, the relative ratio (F y,3D − F y,1D )/F y,3D increases with increasing the ratio E p /E m . In addition, the relative ratio increases with increasing the Poisson ratios ν p and ν m of the particles and matrix, respectively. The relative ratio increases relatively slowly within the interval ν m , ν p ∈ [0.0, 0.2] and sharply when ν m , ν p ≤ 0.38.
The obtained relative ratios may be treated as too big; however, the effective mechanical properties have to know to approximate a particulate composite as a homogeneous solid by the springs. This prognosis is always inaccurate due to many factors affecting the properties of a composite that cannot be taken into account in the calculations. For example, the well known Hashin-Shtrikman lower and upper bounds [14] may also vary within wide intervals.
When Poisson's ratio ν p = ν m = 0, then the obtained axial forces F y,3D of the 3D FE model of the loading case LC1 can be compared with the axial forces F y,HS,low = U y E e f f ,low A/h c and F y,HS,up = U y E e f f ,up A/h c of the homogeneous cube whose effective elastic moduli E e f f ,low and E e f f ,up are calculated by Hashin-Shtrikman's bounds, where A = b c d c is the cross-section area of the cube. The dimensions of the homogeneous cube are the same as for the 3D FE model shown in Figure 6.
The moduli E e f f ,low and E e f f ,up are calculated by taking into account the volume fractions of the particles φ p of the 3DFE model: φ p = 28.46% as d p = 10 mm, φ p = 51.45% as d p = 12 mm, and φ p = 61.35% as d p = 13 mm. The analysis has shown that F y,HS,low is closer to F y,3D than F y,HS,up to F y,3D . It is also determined that the relative ratios |1 − F y,HS,low /F y,3D | < |1 − F y,1D /F y,3D | when φ p ∈ {28.46, 51.45}% or when d p ∈ {10, 12} mm; however, when φ p = 61.35% or d p = 13 mm, then |1 − F y,HS,low /F y,3D | > |1 − F y,1D /F y,3D |. Therefore, when the distance between the particles' surfaces decrease, the efficiency of the proposed methodology increases in comparison to the Hashin-Shtrikman bounds. It should be noted that the already existing methodologies for predicting the effective mechanical properties of composites cannot be applied directly to the lattice model, since the stiffnesses of the connecting element are to be determined. The extra methodology has to be involved in the calculations. These considerations show that the proposed methodology can be useful for predicting the stiffness constants of the connecting element. Moreover, as the analysis showed, the greater the ratio E p /E m is, the more accurate is the proposed methodology in comparison to well-known methodologies of the prediction of the effective mechanical properties of the composites due to the fact that the relative ratios F i,3D and F i,1D , i ∈ {x, y} do not increase very much with increasing the ratio E p /E m .

Conclusions
The evaluation of the axial stiffness of the springs for elastic particulate composites of spherical particles is the main subject of the present article. The methodology takes into account spherical surfaces of the particles and the bond. The closed-form solutions of two upper and lower bounds of stiffness have been obtained.
The obtained formulas have been verified by the three-dimensional FE model of the connecting element. In the majority of cases, the stiffness of the connecting element obtained by the three-dimensional FE model is between the lower and upper bounds. In the analysis, it has also been shown that the lower limit of the stiffness is closer to the results obtained by the three-dimensional FE method than the upper limit when the distance between the particle surfaces is small.
The spring model has been verified by a three-dimensional FE model of the elastic particulate composite cube subject to the tension and the shear force when the lower bound of the axial stiffness of the springs was taken into the calculations. It has been shown in the analysis that the absolute values of the relative ratios of the tensile and shear forces of the spring model and the three-dimensional FE models may reach even up to 42%. However, the prognosis of the effective mechanical properties of the composites using the existing methodologies is almost always inaccurate.. Therefore, the proposed methodology of the evaluation of stiffness of the springs has acceptably high accuracy and can be applied to the spring method of particulate composites. The proposed methodology can also be suitable to evaluate the stiffness of the springs of particulate composites of bonded particles when the diameter of the bond is the same as the diameter of the particles.
Author Contributions: D.Z. participated in: conceptualization, methodology creation, validation of the results, writing the original draft along with editing and reviewing it, and visualization; V.R. participated in: methodology creation, numerical analysis, and validation of the results.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations and main notations are used in this manuscript: D p and D m are elastic constants of the particle and interface member (matrix), respectively E p and E m are moduli of elasticity of the particle and interface member (matrix), respectively F x,1D and F y,1D are shear (in direction x) and normal (in direction y) forces of the the spring model of the particulate composite cube F x,3D and F y,3D are shear (in direction x) and normal (in direction y) forces of the the 3D FE model of the particulate composite cube K s is stiffness of a spring or a connecting element K s,I and K s,I I are limit axial stiffnesses of the spring K s,FEM is stiffness of the connecting element obtained by the FE model LC1, LC2 are loading cases L c is distance between particles' surfaces SM is a spring model consisting of pin-connected springs U x and U y are displacements in directions x and y, respectively R p and d p are radius and diameter of a particle R b is radius of the interface member ν p and ν m are Poisson's ratios of the particle and interface member (matrix), respectively φ p is particle volume fraction of a composite