Analysis of Boundary Layer Influence on Effective Shear Modulus of 3-1 Longitudinally Porous Elastic Solid

Corresponding author https://doi.org/10.1590/1679-78255960 Abstract The evaluation of the effective properties of nonhomogeneous solids using analytical methods is, in general, based on the assumption that these solids have infinite dimensions. Here, we investigate the influence of both the number of holes and the boundary layer of a solid with finite dimensions on the determination of these properties. We use the Asymptotic Homogenization Method (AHM) to determine the effective shear modulus of an elastic solid with infinite dimensions containing a uniform and periodic distribution of circular cylindrical holes arranged on a hexagonal lattice. We also use the Finite Element Method (FEM) to determine this modulus in the case of a solid with finite dimensions containing the same uniform distribution of cylindrical holes away from its boundary. Near the boundary, we consider a layer of material with no holes, which is usually left in the fabrication process of samples. Both solids have the same elastic properties and are subjected to similar anti-plane shear loadings. For the finite medium, we study two sequences of domains discretized by the FEM, which are called the Fixed Layer Sequence (FLS) and the Fixed Domain Sequence (FDS). For the FLS, the layer thickness is kept fixed and both the dimensions of the domain and the number of holes vary. For the FDS, the dimensions of the domain are kept fixed and both the number of holes and the layer thickness vary. Results obtained from numerical simulations are then used to generate graphs of the effective shear modulus versus void volume fraction. It is observed that, in the FLS case, the shear modulus obtained from the numerical simulations converges to the analytical solution obtained via AHM. It is also observed that, in the FDS case, the shear modulus obtained from the numerical simulations converges to a limit function, which is close to the analytical solution obtained via AHM. For comparison purposes, we have also calculated the effective shear modulus of porous elastic solids containing a square array of circular cylindrical holes. We then show graphs of this modulus versus void volume fraction for both hexagonal and square arrangements that are very close to each other up to void volume


INTRODUCTION
The study of the behavior of nonhomogeneous solids, such as biological tissues, requires the determination of their effective properties via analytical methods, computational methods, or a combination of both. Generally, in the case of analytical methods, we assume that the solid has infinite dimensions. Since the samples used in laboratory experiments have finite dimensions, it is of interest to verify whether such expressions provide accurate values for the constants obtained via these experiments. A step towards this goal is to simulate these experiments numerically.
Several numerical approaches have been proposed to predict the overall constitutive behavior of nonhomogeneous solids, among which we cite Ramakrishnan and Arunachalam (1990), Said et al. (2016), Pingaro et al. (2019), Ghosh et al. (1996), Ghosh et al. (2001), Marino et al. (2019), andLo Cascio et al. (2020), who have considered random distributions of heterogeneities. These authors study how spatial distributions of heterogeneities and volume fraction influence the effective elastic properties of nonhomogeneous solids by identifying a Representative Volume Element (RVE) and determining the proper boundary conditions that should be applied on the boundary of the RVE. Different from these works, here, we apply conditions on the boundary of the whole domain.
Some works related to the evaluation of effective properties in cortical bones are discussed below. Swan et al. (2003) neglect the lamellar structure of cortical bone matrix and model the resulting Haversian bone matrix as a homogeneous and isotropic linearly elastic medium. They consider a unit cell model of Haversian bone with squarepacked non-overlapping osteons. The cell has either or Haversian porosity and is submitted to periodic boundary conditions. They use the Finite Element Method (FEM) to compute all the effective poroelastic moduli of this cell. In this work, we also neglect the lamellar structure of the cortical bone matrix but, differently from these authors, we use a hexagonal array of cells containing unidirectional circular cylindrical holes to model non-overlapping osteons and we employ finite element discretization of the whole domain, instead of the unit cell model, to compute the shear effective modulus. Grimal et al. (2008) introduce a method to obtain the elastic tensor of a cortical bone sample at the mesoscale from a mapping of its microscale elasticity. The mesoscale properties are estimated based on a finite element homogenization procedure and the results are compared with available experimental data. The authors observe, however, that, although experimental data indicate that always hold, their computed effective elastic coefficients do not obey these expressions, prompting them to question the validity of their computations for the case of effective We have verified that by using both analytical and computational methods. Parnell and Grimal (2009) use the asymptotic homogenization method (AHM) to predict the influence of porosity on the induced anisotropy of cortical bone. For this, they assume that the matrix phase of the bone is made of a homogeneous and isotropic linearly elastic material. For the case of circular pores, the authors show that the results obtained via AHM are in good quantitative and qualitative agreement with finite element results from simulations of problems using real two-dimensional microstructures obtained from images of cortical bone. Later, Parnell et al. (2012) compare theoretical predictions of the effective elastic moduli of cortical bone at both the meso-and macroscales. They consider the efficacy of three alternative homogenization approaches: the AHM, the Mori-Tanaka scheme and the Hashin-Rosen bounds. The authors point out that, although the mesoscale behavior of bone is widely accepted as important, models incorporating its effect have started to appear only recently. Except for the work of Swan et al. (2003), the authors above working on cortical bones have considered circular cylindrical holes arranged on a hexagonal lattice.
In addition to the above works on homogenization schemes to evaluate effective properties of cortical bone, the impact of the elastic symmetry, assumed for the bone matrix, on the mesoscopic behavior was discussed by Sevostianov and Kachanov (2000), among others. These authors have proposed a micromechanics model to study the influence of porosity on the anisotropy of cortical bone. They have concluded that the differences between the cases of empty pores and pores filled with soft material were insignificant. Grimal et al. (2008) have quantified this influence and concluded that, on the contrary, it is significant. In the face of these disagreements and lack of methodologies in the literature to associate the complex interaction between bones' structural and mechanical properties, we decided to model only solids containing empty holes on the application of our numerical and analytical results.
The aim of this research is then to develop a reliable and computationally efficient method of analysis suitable for predicting the response of elastic solids containing empty holes and the influence of their boundaries on the determination of their material properties. Specifically, we employ the Asymptotic Homogenization Method (AHM) and the Finite Element Method (FEM) to determine the effective properties of an elastic solid that has a uniform and periodic distribution of circular cylindrical holes in an isotropic linearly elastic medium. The cylinders are centered in unit cells of hexagonal cross sections. In order to apply the AHM, we consider that the microstructure of the solid consists of two phases distributed periodically over a domain that has infinite dimensions. A cross section of the solid is shown in Fig. 1 This distribution of phases allows us to expand the solution of the related equilibrium problem in terms of an asymptotic series and obtain local problems. The solutions of these local problems can be calculated analytically for solids with simple microstructures, such as the one considered in this work, or, more generally, numerically to include solids with complex microstructures, such as bone. These local solutions are then used in the calculation of the effective modulus of the elastic solid, which depends upon the physical and geometrical properties of its phases.
Specifically, we use closed form expressions obtained by Bravo-Castillero et al. (2009) and Aguiar et al. (2013) to calculate analytically the effective shear elastic modulus of the solid. The procedure used in this calculation can be used to evaluate the other effective elastic moduli of the solid. The procedure may also be used to evaluate the elastic moduli of a solid with different microstructures, such as the one with unit cells having square cross sections, which is also used in this work for comparison purposes. Other more general approaches can also be found in the literature. In particular, López-Realpozo et al., (2011), Otero et al., (2013, López-Realpozo et al. (2014), and Otero J. A. et al. (2016) have used the AHM to formulate local problems on unit cells with periodic boundary conditions and solved these problems either analytically or computationally using FEM. These solutions were used to calculate the effective properties of two-phase fibrous periodic composites with imperfect contact conditions at the interface for simple geometries, such as circular fibers distributed in a rhombic array. This approach is different from ours, mainly because we formulate our problem on the whole domain, and not on a unit cell.
Next, we simulate numerically laboratory experiments using FEM and we use the results of these simulations to evaluate the effective elastic modulus . For this, we consider two sequences of computational domains, which are called the Fixed Layer Sequence (FLS) and the Fixed Domain Sequence (FDS). For the FLS, the layer thickness is kept fixed and both the dimensions of the domain and the number of holes vary. For the FDS, the dimensions of the domain are kept fixed and both the number of holes and the layer thickness varies. To the best of our knowledge, this is the first study using the FDS. We then compare the numerical evaluation of with the corresponding analytical evaluation obtained through AHM. Recall from above that the solid has infinite dimensions in the AHM analysis and finite dimensions in the numerical simulation. The numerical experiment consists of the anti-plane shear of a cylindrical sample with rectangular cross section containing the same uniform distribution of unit cells of the infinite solid away from the boundary of the solid. Near the boundary, we consider a layer of material and investigate the influence of this layer on the evaluation of the effective modulus. The cross section of the solid is illustrated in Fig. 2. The material of this sample is the same isotropic linearly elastic material considered in the analytical approach.
This investigation represents an ongoing effort of the research group to obtain the effective moduli of elastic solids analytically and numerically, and it will be useful to evaluate the influence of the boundary layer in 3D printed samples for laboratory testing. One of the most important 3D printing techniques is Additive Manufacturing (AM). Recently developed, AM is an emerging technological advancement in the field of engineering. See, for instance, Yap et al. (2019) and Wohlers (2012), due to its highly adaptable manufacturing capabilities and ease of use, for instance, for rapid prototyping of intricate geometries with metals, polymers, and fiber-reinforced composite materials. In AM, the sample boundary layer is comprised of several shell perimeters, which are among the main control factors on the mechanical properties of samples fabricated in, for instance, polylactic acid (PLA). According to Lanzotti et al. (2015), other important control factors are the layer thickness of deposition and the infill orientation of each layer. The high variability of results observed by these authors together with a shortage of existing literature concerning the impact of these factors on the mechanical properties of samples made by additive manufacturing has motivated this work. They investigated the main impact of three process parameters -layer thickness, infill orientation and the number of shell perimeters -on the mechanical properties of parts fabricated in PLA and the effects of interactions. They use a response surface methodology to propose an empirical model, which connects process parameters and mechanical properties.
In another work, Ćwikła et al. (2017) study the selected mechanical strength properties of 3D printed elements carried out on a set of standardized samples. In their work, Ćwikła et al. (2017) show that the shell thickness has a key influence on the tensile strength of the samples; for example, if the maximum strength is the priority, shell thickness should be increased. These authors observe that for increasing values of shell thickness over , the infill of relatively small samples is practically replaced by the solid, closely extruded filament threads. These threads have a much higher tensile strength than the standard infill pattern, resulting in overall tensile strength similar to infill specimens. Ćwikła et al. (2017) mention also that the results indicate the need for further research concerning, e.g., simultaneous changes in many parameters and rigidity of obtained elements, which requires the development of specimens appropriate for the specific study of 3D printed objects, as well as methods of measurement. were considered: material, infill density, infill pattern, and outer vertical shell. The authors test three FFF distinct materials (polycarbonate (PC), acrylonitrile butadiene styrene (ABS), and polylactic acid (PLA)) and they investigate three infill densities ( , , ) along with two different infill patterns (rectangular and honeycomb). These authors also investigate printing accuracy, repeatability and mechanical behavior of porous 3D printed structures. Results obtained via FEM by them indicate that both PC and ABS can be safely used to fabricate a porous lumbar cage with a honeycomb infill density and a honeycomb infill pattern, concluding that the 3D printing assisted FEM can be used to predict the performance design with varying manufacturing parameters and potentially reduce product design and development time. In this respect, our work contributes to explain how shell thickness influences the effective elastic modulus and will be useful in experimental applications, such as, the one presented in Provaggi et al. (2019).
The paper is organized as follows. In Section 2 we present the statement of the anti-plane shear problem for an isotropic linearly elastic medium containing uniform and periodic distributions of circular cylindrical holes with infinite dimensions arranged on a hexagonal array. In Section 3 we apply this formulation together with the AHM in the calculation of the effective elastic constant . In Section 4 we consider finite dimensions for the elastic medium and we use the analogy between the associated elastic problem and the linear steady state heat conduction problem to formulate the problem that we solve numerically by using a finite element commercial package. In Section 5 we then compare the results obtained via AHM and via FEM by taking into account the number of holes, concentration of holes and the boundary layer thickness of the solid with finite dimensions. In Section 6 we present the conclusions of this work.

PROBLEM STATEMENT
Consider an isotropic linearly elastic medium with infinite dimensions containing a uniform and periodic distribution of circular cylindrical holes centered in unit cells with hexagonal cross sections. With respect to a Cartesian coordinate system with origin at the axes of the holes are parallel to the coordinate axis and the cross section of the region occupied by the medium in a stress-free undeformed configuration is parallel to the plane . This cross section together with the cross section of the hexagonal unit cell are shown in Fig. 1, where and correspond to the solid portions of the medium, without holes, and and correspond, respectively, to the radius and the boundary of the hole in the cell. The medium is in equilibrium in the absence of body force and it is subject to loading that is orthogonal to the plane at infinity, yielding a linear shear antiplane problem. In this problem, we use Voigt notation for the components of the elasticity tensor, in which we employ the usual simplification on the notation of the moduli taking the subscript indices to . In this work, the components of interest are and .
We consider that the only nonzero component of the displacement field is parallel to the -direction and depends on the coordinates only, i.e., . Substituting these components in the straindisplacement relations (1) and the resulting expressions into the generalized Hooke's law for isotropic materials, we obtain that the only nonzero components of the stress tensor are given by (2) in which we recall from above that are the elastic moduli. On the other hand, the only nonzero equilibrium equation with no body force is given by (3) Substituting the constitutive relations in Eq. (2) into the equilibrium equation given by Eq. (3) yields the differential governing equation (4) where is the two-dimensional Laplace operator with respect to the global variable Next, the stress-free condition on the walls of the holes, which we denote by , is given by (5) where are the components of the normal to a cylindrical surface.
Substituting the constitutive relations in Eq.
(2) into the boundary conditions in Eq. (5), we obtain (6) where is the gradient operator with respect to the variable The first boundary-value problem of interest in this work consists of finding the displacement field that satisfies the differential equation in Eq. (4), the boundary condition in Eq. (6), and periodic conditions on at infinity.

THE HOMOGENIZATION PROCEDURE
In this section, we review basic aspects of the AHM, which are used to obtain a closed-form expression for the effective elastic modulus of the elastic medium with infinite dimensions. Further details can be found in Aguiar et al. (2013). This expression is used in Section 5 for comparison purposes with numerical approximations of obtained for a medium with finite dimensions.
Let be a small geometric parameter defined as the ratio between a characteristic length of the cell, such as the length of a side of the cell, and a characteristic length of the region occupied by the two-phase composite.
For a refined microstructure, and we employ a two-scale expansion method to expand the displacement field as a power series of in the form 1 , , 1, 2, 3, 2 Latin American Journal of Solids and Structures, 2020, 17(8 Thematic Section), e313 6/19 (7) in which is a local variable depicted in Fig. 1.b, as opposed to the global variable , are twicedifferentiable functions of and represents the displacement of a homogenized body. Moreover, are twice-differentiable functions of both and are highly oscillating, and represent correction terms to the zerothorder approximation .
(2) and noting that is a function of the local variable which vanishes inside the holes and is constant outside, we obtain (8) where (9) . Substituting the expressions in Eq. (8) together with Eq. (9) into Eq. (3), adopting the Einstein summation convention, and applying the chain rule yields (10)   (11) The governing equations in the homogenized medium can be obtained by employing the averaging operation to the equations Eq. (11) and using the periodicity of in Thus, it yields (12) in which Next, we substitute the equations Eq. (9) into Eq.(10), yielding Recalling from Eq. (7) that depends on only the variable and integrating both sides of the Eq. (13) with respect to , we proceed using the separation of variables taking (14) where are twice-differentiable functions of and satisfy periodic conditions on the external boundary of Substituting the expressions in Eq. (14) into Eq. (13) (15) where (16) The differential equations in Eq. (15) are used to determinate the Y-periodic functions , below.
Substituting the expressions in Eq. (14) into Eq. (9), we obtain (17) Applying the average operator, defined below Eq. (11), onto Eq. (17), we obtain (18) where (19) is the general expression of the effective elastic constant of the homogenized medium. Now, to calculate the effective properties using Eq.  (20) where is the two-dimensional Laplace operator with respect to the local variable Substituting Eq. (8) in Eq. (5) yields and letting we obtain the zeroth-order term . On the other hand, substituting expansions (7) into the constitutive relations (2) and substituting the resulting expressions together with the general forms (14) into the boundary conditions (5)   For the problem at hand, which is to calculate analytically the effective shear elastic modulus in the longitudinal direction, of the elastic solid, only the Local Problem over the local domain, illustrated on the right-hand side of Fig. 1, will be specified at the outset and solved. This problem consists of finding a function from Eq. (20)-(22), harmonic and of zero average in that satisfies the system of equations (23)   (24)   (25) where is the two-dimensional Laplacian in corresponds to the boundary of the hole in the cell, and with being the volume of . Nevertheless, by solving Local Problem , we verify that the effective shear elastic modulus is equal to .
The solution of Eq. (23)-Eq. (25) is sought in the class of doubly periodic harmonic functions, , which depends on the complex variable , in the form of a series with undetermined real coefficients , given by (26) where is the angle of the unit cell, which is for the hexagonal cell, and and are, respectively, the quasi-periodic Weierstrass zeta function and its kth derivatives of periods and . Observe from Eq. (23)-Eq. (25) that the local solution is not dependent of the physical properties of the material; it depends on the geometry and the arrangement of the holes.
Once the local problem is solved, the effective modulus is given by (27) where is the shear modulus of the solid part of the elastic medium and it can be shown (Bravo-Castillero et al., 2009) that (28) where is the solution of the infinite system (29) where is the Kronecker delta and ef 44 , c

COMPUTATIONAL METHOD
Let be a rectangular region containing a part of the periodic distribution of hexagonal unit cells considered in Section 2 and let be a thin layer of thickness surrounding , as illustrated by the blue frame in Fig. 2. The layer has the same material properties of the solid part occupying . The whole region has dimensions . Let also be the whole region excluding the holes, in analogy to the region of Section 2. The region is occupied by the same homogeneous, isotropic, and linearly elastic material considered in Section 3. On the external boundary of the sample, we impose zero displacement on the left-hand side, , zero tractions on the lower, , and upper, , sides, and displacement on the right-hand side, , which has magnitude and is applied normal to the -plane. Also, is the outward normal vector to the external boundary.
12 12 (1 6 1) (6 1 6 1) (6 1 6 1) ( 6 1 6 1)  1 , , To solve numerically the equilibrium problem described at the end of in Section 2 in the domain , we observe that this problem is analogous to a linear steady state heat conduction problem, which consists of finding the temperature field that satisfies where and are, respectively, the Laplacian and the gradient operators defined in , is the union of all the contours of the holes, and is the external boundary of . We use the finite element commercial package COMSOL 4.4 ® to obtain approximate solutions of Eq. (33)-Eq. (35). The temperature field is approximated by secondorder triangular Lagrange elements using free mesh generation automatic triangulation, and a direct finite element solver.
Once an approximate solution of the thermal problem is obtained for a given discretization, we associate the temperature with the displacement and the conductivity with the elastic modulus . The effective elastic modulus can then be determined from the expression (36) where is the heat flux on .
In order to investigate the influence of both the dimensions of the domain and the number of holes on the convergence of the numerical solutions obtained via FEM to analytical solutions obtained via AHM, we create two sequences of discretized domains. In the first sequence we vary the dimensions of the domain and the number of holes and keep fixed the thickness of the layer represented by the blue frame in Fig. 2. This approach we call Fixed Layer Sequence (FLS). This approach is based on a similar approach used by Hollister and Kikuchi (1992), who have fixed the dimensions of a cell, understood as the basic repeating unit of a periodic material structure, and piled them up to build samples. The authors compare results obtained from a homogenization theory to results obtained from a finite element implementation of standard mechanics of materials approaches for the analysis of periodic porous composites. The authors analyze two-dimensional cellular structures with solid volume fractions of 30%, 50%, 70%, and 90%. For the standard mechanics approach, RVEs containing 1, 4, 9, 16, and 25 cells were analyzed for each solid volume fraction. In the second sequence we keep the dimensions of the domain fixed and vary both the number of holes and the thickness . We call this approach Fixed Domain Sequence (FDS). We have not found a similar approach in the literature.
In Figs. 3 and 4 we show computational domains belonging to FLS and FDS, respectively. To construct Fig. 3, first, observe from the right-hand side of Fig. 2 that the unit cell of hexagonal shape having area of 10 -4 m 2 yields h = 0.0062 m. We then consider both a volume fraction of the holes and a layer thickness a = h/2 = 0.0031 m. A sample is then constructed by repeating this unit cell. In this work we vary the dimensions of in Fig. 2 according to the cases (a) , (b) , (c) , (d) . In Fig. 4 the dimensions of are fixed and given by , the layer thickness is , the volume fraction of the holes is , and we vary the dimension of the hexagon according to the cases (a) h = L/4, (b) h = L/8, (c) h = L/16, (d) h = L/32. Because of the way the samples were constructed, which reflects limitations experienced during the generation of the geometry, the computational domain has a rectangular shape in the FLS approach and a square shape in the FDS approach. Below, we show computational results indicating that the shape of the geometry of has no influence on the convergence analysis.  In Section 5 we present results for the effective modulus obtained from these discretizations. In the first two cases, corresponding to the values , we compare our results with results obtained by Swan et al. (2003). For these two cases, observe from Fig. 5 that the walls of the holes do not intersect the boundary layer .
Intersections occur for near . The case of is then introduced as an intermediate case to show these intersections. Finally, represents the limit case where the walls of the holes are nearly touching each other.

RESULTS AND DISCUSSION
Here, we present and discuss some analytical and computational results concerning the determination of the effective elastic modulus . Specifically, we use the FDS and FLS approaches discussed in Section 4 to investigate convergence of computational results obtained via FEM to analytical results obtained via AHM. For this, we consider as the elastic modulus of the solid occupying either the region , which has infinite dimensions, or , which has finite dimensions.
First, we use the FDS approach, with fixed dimensions of given by and fixed layer thickness given by , and consider a sequence of decreasing values of the hexagon dimension in the set Except for h = L/2, the other values of correspond to the cross sections in Fig. 4. In Fig. 7 we show graphs of the effective modulus versus the void volume for this sequence. In the figures below, the label stands for a Finite element result obtained with a layer thickness having +1 significant digits and using the FDS approach. Thus, F3 # in Fig. 7 means that we have used the FDS approach together with . Recall from Section 4 that these graphs were obtained via FEM together with Eq. (36). We also show a curve obtained via AHM by using Eq. (31) together with Eq. (32). Note that the curves obtained via FEM converge non monotonically to the curve obtained by AHM as , even at a high void volume fraction . These results show the fact that the numerical result for this constant is influenced by the decreasing values of the hexagon side dimension in the samples. This convergence is expected from the AHM theory summarized in Section 3, according to which the displacement field is expanded in power series of the small parameter in Eq. (7), which, here, can be taken as either or .

Figure 7
Effective elastic modulus versus void volume fraction using the FDS approach depicted in Fig. 4.
Next, we compare the best numerical curve presented in the Fig. 7, corresponding to , with the curve obtained via AHM and the best curve by using FLS approach, corresponding to obtained by Aguiar et al. (2018). In the figures below, the label , without the superscript #, stands for a Finite element result obtained with the FLS approach together with the ith cross section in Fig. 3. Thus, F4 means that the cross section (d) in Fig. 3 was used. For this case, the authors considered elastic samples with rectangular cross sections having the fixed dimensions and the thickness of the boundary layer . This figure shows that for a thin boundary layer fixed we obtain good agreement with the analytical solution, obtained via AHM, for two independent methods.  Thus, in Fig. 8 we show again the curves obtained via AHM and via FEM with . These curves were obtained numerically via FEM by using two independent methods, the FDS and FLS approaches, and compared with that one obtained analytically by AHM. Observe from this figure that all curves are very close to each other. This result allows us to choose as a representative parameter of number of holes in a fixed domain and below we will vary the boundary layer thickness Next, we compare results obtained from both approaches FDS and FLS. In the FDS approach, we still take , but hold the side of the hexagon fixed at and consider a sequence of decreasing values of the layer thickness in the set . For the FLS approach, we consider the data and the corresponding discretizations shown in Fig. 3. In Fig. 9 we show graphs of the effective modulus versus the void volume for both sequences. In Fig. 9(a) we consider the whole interval and in Fig. 9(b) we show a zoom in of the graphs for the interval . Again, these graphs were obtained via FEM together with Eq. (36). Also, the graphs obtained from the FLS approach were originally obtained by Aguiar et al. (2018) and are reproduced here for comparison purposes. Observe from Fig. 9(a) that the numerical results obtained from the FDS approach, which are represented by , converge to limit values that are close to the analytical results obtained from AHM as decreases. In the case of the FLS approach, there is convergence to the analytical results, which means that the bigger the rectangular cross sectional area the closer the numerical solution is to the analytical solution. This convergence behavior is also expected from the AHM theory, because the parameter in Eq. (7), which, again, can be taken as either or , decreases as both and increase for a fixed . In particular, good convergence is observed at all values of void volume fraction shown in Fig. 9 (b). In addition, near the limit point where the walls of the holes touch each other, which we recall from Section 4 is approximately given by and corresponds to the cases (c) and (d) of Fig. 6, no convergence was obtained.
In Fig. 10     In addition to the calculation of the effective modulus for solids containing a hexagonal array of holes, we also calculated the effective properties of solids containing a square array of holes. In Fig. 11 we show graphs of versus obtained from AHM and both FLS and FDS approaches using square (SQ) and hexagonal (HEX) arrays of holes. The labels SQ AHM and HEX AHM correspond to AHM results obtained with square and hexagonal arrays of holes, respectively. Analogous correspondence holds between the labels containing F4 and F4 # and the FLS and FDS results, respectively. Fig. 11(a) refers to the whole interval (0,1) and Fig. 11(b) refers to the interval (0,0.04). These figures show that all the curves are indistinguishable for volume fractions up to 0.5. We then see from these figures that both arrays of holes yield the same values of . Next, we compare our results with the results obtained by Swan et al. (2003) for the corresponding effective shear elastic constant using the same properties given by these authors. In that work, however, the holes are centered in a square unit cell (instead of the hexagonal cell). In Table 1 we present the void volume fraction , which has the values and , the effective shear elastic coefficient obtained by Swan et al. (2003) from the fully drained case, the  Table 2 together with their Table 1 that is associated with Young's modulus and Poisson's ratio and that is associated with and . The fully drained case corresponds to the response of bone in which the fluid carries no excess pressures from applied loadings. Observe from Table 1 a very good agreement between our calculations and theirs for the case and a larger discrepancy in the case , which is being investigated. In this paper we have introduced the FDS method, which is a reliable and computationally efficient method to build samples for both the numerical prediction of the response of elastic solids containing empty holes and the evaluation of the influence of the boundaries on the determination of their material properties. The research was focused on the calculation of the effective modulus in linear elastic solids with two scales only. Our numerical method may, however, be further developed to consider multiscale and hierarchical structures similar to the hierarchical structure studied by Ramírez-Torres et al (2018). These authors investigate the effective properties of hierarchical composites at each structural level and apply the results of this investigation in the analysis of linear elastic composites with hierarchical structure, where the calculated effective properties at the lower structural level become the known elastic properties for the problems arising at the higher level. In addition to the investigation of multiscale and hierarchical structures, we are also interested in using ideas based on the FLS and FDS methods to evaluate the influence of the boundary layer in 3D printed samples for laboratory testing. The goal is to compare the effective constants of these porous solids with the corresponding constants obtained analytically via AHM and numerically via FEM.

CONCLUSION
In this paper, we have shown some results concerning the analytical and numerical evaluations of the effective modulus of linear elastic solids having a periodic distribution of voids arranged in hexagonal cells. When the side of the hexagon with length decreases and the boundary layer thickness is fixed, or, when decreases and is fixed, we observe convergence of obtained computationally via Eq. (36) to the corresponding modulus obtained analytically via both Eq. (31) and Eq. (32). Another result of interest is that both arrangements, hexagonal and square, yield effective elastic modulus that are very close to each other up to for both and tending to zero. This investigation represents an ongoing effort to obtain the effective moduli of elastic solids using analytical and computational methods. Future work includes obtaining these moduli experimentally so that we can compare them with both the analytical and computational predictions.