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

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 fraction of 0.5.

Keywords:
linear elasticity; asymptotic homogenization method; finite element method; effective modulus; boundary layer

Graphical Abstract

1 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), and Lo 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 square-packed non-overlapping osteons. The cell has either $1\%$ or $4\%$ 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 $c_{66}<c_{55}=c_{44}$ 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$c_{44}$ We have verified that $c_{66}^{ef}\le c_{55}^{ef}=c_{44}^{ef}$ 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 $c_{44}^{ef}$. 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 $c_{44}^{ef}$ 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 $4$, 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 $100\%$ 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.

Provaggi et al. (2019) investigate the use of Fused Filament Fabrication (FFF) together with FEM to the understanding of manufacturing parameters in the design process of a lumbar fusion cage. The following parameters 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 ($25\%$, $50\%$, $75\%$) 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 $50\%$ 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 $c_{44}^{ef}$ 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$c_{44}^{ef}$. 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.

2 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 $O,$ the axes of the holes are parallel to the coordinate axis $O{x}_3$ and the cross section of the region occupied by the medium in a stress-free undeformed configuration is parallel to the plane $x_1{x}_2$. This cross section together with the cross section of the hexagonal unit cell are shown in Fig. 1, where $M$ and $Y$ correspond to the solid portions of the medium, without holes, and $R$ and $\Gamma$ correspond, respectively, to the radius and the boundary of the hole in the cell.

Figure 1
Cross section of a 3-1 longitudinally porous elastic medium of infinite dimensions (left) and the corresponding periodic cell for a hexagonal array. Adapted from Bravo-Castillero et al. (2009).

The medium is in equilibrium in the absence of body force and it is subject to loading that is orthogonal to the plane $x_1{x}_2$ 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$c_{ijmn},$ taking the subscript indices $ijmn,$ $i,j,m,n=\mathrm{1,2,3,}$ to $11\to \mathrm{1,}22\to \mathrm{2,}33\to \mathrm{3,}23\to \mathrm{4,}13\to 5and12\to 6$. In this work, the components of interest are $c_{1313}=c_{55}$ and $c_{2323}=c_{44}$.

We consider that the only nonzero component of the displacement field $u:M\to {ℝ}^3$is parallel to the $x_3$-direction and depends on the coordinates $(x_1,x_2)$ only, i.e., $u:=u_3^{}(x_1,x_2)$. Substituting these components in the strain-displacement relations

${\epsilon }_{ij}^{}=\frac{1}{2}\left(\frac{\partial u_i^{}}{\partial x_j}+\frac{\partial u_j^{}}{\partial x_i}\right),i,j=\mathrm{1,2,3,}$ (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

${\sigma }_{i3}^{}=c_{44}^{}\frac{\partial u}{\partial x_i},i=\mathrm{1,2,}$ (2)

in which we recall from above that $c_{44}^{}=c_{55}^{}$ are the elastic moduli. On the other hand, the only nonzero equilibrium equation with no body force is given by

${\sigma }_{i3}^{}=c_{44}^{}\frac{\partial u}{\partial x_i},i=\mathrm{1,2,}$ (3)

Substituting the constitutive relations in Eq. (2) into the equilibrium equation given by Eq. (3) yields the differential governing equation

${\Delta }_{x}u:=\frac{{\partial }^{2}u}{\partial x_1^2}+\frac{{\partial }^{2}u}{\partial x_2^2}=0inM,$ (4)

where ${\Delta }_x:=\left({\partial }^2/x_1^2\right)+\left({\partial }^2/x_2^2\right)$ is the two-dimensional Laplace operator with respect to the global variable $x.$

Next, the stress-free condition on the walls of the holes, which we denote by $\partial M$, is given by

${\sigma }_{i3}^{}{n}_i=\mathrm{0,}i=\mathrm{1,2,}$ (5)

where $n_i$are the components of the normal $n$ to a cylindrical surface.

Substituting the constitutive relations in Eq. (2) into the boundary conditions in Eq. (5), we obtain

$n\cdot {\nabla }_{x}u:=n_1\frac{\partial u}{\partial x_1}+n_2\frac{\partial u}{\partial x_2}=0on\partial M,$ (6)

where ${\nabla }_x:=\left(\partial /x_1^{}\right)e_1+\left(\partial /x_2^{}\right)e_2$ is the gradient operator with respect to the variable $x.$

The first boundary-value problem of interest in this work consists of finding the displacement field $u:M\to ℝ$ that satisfies the differential equation in Eq. (4), the boundary condition in Eq. (6), and periodic conditions on $M$ at infinity.

3 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 $c_{44}^{ef}$ 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 $c_{44}^{ef}$obtained for a medium with finite dimensions.

Let $\epsilon$ 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, $\epsilon <<\mathrm{1,}$ and we employ a two-scale expansion method to expand the displacement field as a power series of $\epsilon$ in the form

$u_{}^{\epsilon }(x)=u_{}^{(0)}\left(x\right)+\epsilon u_{}^{(1)}\left(x,y\right)+{\epsilon }^2{u}_{}^{(2)}\left(x,y\right)+\cdots ,$ (7)

in which $y=\left(y_1,y_2\right)$ is a local variable depicted in Fig. 1.b, as opposed to the global variable $x$, $u_{}^{(n)},n=\mathrm{0,1,2...}$ are twice-differentiable functions of $x,$ and $u_{}^{(0)}$represents the displacement of a homogenized body. Moreover, $u_{}^{(n)},n=\mathrm{1,2,3,...}$are twice-differentiable functions of both $x$ and $y,$ are highly oscillating, and represent correction terms to the zeroth-order approximation $u_{}^{(0)}$.

Substituting Eq. (7) into Eq. (2) and noting that $c_{44}^{}=c_{55}^{}$ is a function of the local variable $y,$ which vanishes inside the holes and is constant outside, we obtain

${\sigma }_{i3}\left(x\right)={\sigma }_{i3}^{(0)}(x,y)+\epsilon {\sigma }_{i3}^{(1)}(x,y)+{\epsilon }^2{\sigma }_{i3}^{(2)}(x,y)\dots ,$ (8)

$i=\mathrm{1,2,}$ where

${\sigma }_{i3}^{(n)}(x,y)=c_{44}(y)\left[\frac{\partial u_{}^{(n)}(x,y)}{\partial x_i}+\frac{\partial u_{}^{(n+1)}(x,y)}{\partial x_i}\right],$ (9)

$n=\mathrm{0,1,2,3,...}$. Substituting the expressions in Eq. (8) together with Eq. (9) into Eq. (3), adopting the Einstein summation convention, and applying the chain rule yields

$\frac{\partial {\sigma }_{i3}^{(0)}(x,y)}{\partial y_i}=\mathrm{0,}$ (10)

$\frac{\partial {\sigma }_{i3}^{(0)}(x,y)}{\partial x_i}+\frac{\partial {\sigma }_{i3}^{(1)}(x,y)}{\partial y_i}=0.$ (11)

The governing equations in the homogenized medium can be obtained by employing the averaging operation $\langle \mathrm{...}\rangle \triangleq {\displaystyle \underset{Y}{\int }\left(\mathrm{...}\right)dy/\left|Y\right|}$ to the equations Eq. (11) and using the periodicity of ${\sigma }_{ij}^{(1)}\left(x,y\right)$ in $y.$ Thus, it yields

$\frac{\partial {\overline{\sigma }}_{i3}^{(0)}\left(x\right)}{\partial x_i}=\mathrm{0,}$ (12)

in which ${\overline{\sigma }}_{i3}^{(0)}\left(x\right)=\langle {\sigma }_{i3}^{(0)}\left(x,y\right)\rangle .$ Next, we substitute the equations Eq. (9) into Eq.(10), yielding

$\frac{\partial }{\partial y_i}\left[c_{44}(y)\frac{\partial u_{}^{(1)}(x,y)}{\partial y_i}\right]=-\frac{\partial c_{44}(y)}{\partial y_i}\frac{\partial u_{}^{(0)}(x,y)}{\partial x_i},i=\mathrm{1,2,}(nosumoveri).$ (13)

Recalling from Eq. (7) that $u^{(0)}$ depends on only the variable $x,$ $u_{}^{(0)}(x,y)=u_{}^{(0)}(x),$ and integrating both sides of the Eq. (13) with respect to $y_i$, we proceed using the separation of variables taking

$u_{}^{(1)}(x,y):={}_{l3}U\left(y\right)\frac{\partial u_{}^{(0)}(x)}{\partial x_l},$ (14)

where ${}_{l3}U\left(y\right),l=\mathrm{1,2,}$ are twice-differentiable functions of $y$ and satisfy periodic conditions on the external boundary of $Y.$ Substituting the expressions in Eq. (14) into Eq. (13) yields

$\frac{\partial {}_{k3}\sigma {}_{i3}(y)}{\partial y_i}=-\frac{\partial c_{44}(y)}{\partial y_k},i,k=\mathrm{1,2,}(nosumoverk),$ (15)

where

${}_{k3}\sigma {}_{i3}\left(y\right):=c_{44}(y)\frac{\partial {}_{k3}U{}_{}\left(y\right)}{\partial y_i},i,k=\mathrm{1,2}.$ (16)

The differential equations in Eq. (15) are used to determinate the Y-periodic functions ${}_{k3}U\left(y\right),k=\mathrm{1,2}$, below.

Substituting the expressions in Eq. (14) into Eq. (9), we obtain

${\sigma }_{i3}^{(0)}(x,y)=\left[c_{44}(y)+{}_{k3}\sigma {}_{i3}(y)\right]\frac{\partial u_{}^{(0)}(x)}{\partial x_k},i=\mathrm{1,2}(nosumoveri).$ (17)

Applying the average operator, defined below Eq. (11), onto Eq. (17), we obtain

${\overline{\sigma }}_{i3}^{}(x):=\langle {\sigma }_{i3}^{(0)}\left(x,y\right)\rangle =c_{44}^{ef}\frac{\partial u_{}^{(0)}(x)}{\partial x_i},$ (18)

where

$c_{44}^{ef}=\langle c_{44}(y)+{}_{i3}\sigma {}_{i3}(y)\rangle ,i=\mathrm{1,2}(nosumoveri)$ (19)

is the general expression of the effective elastic constant of the homogenized medium.

Now, to calculate the effective properties using Eq. (19), we need to find the functions ${}_{kl}U{}_{}\left(y\right),$introduced in Eq. (14). These functions are solutions of the equations Eq. (15)-(16). Using the expressions Eq. (15)-(16) and Eq. (17), we obtain

$\Delta {}_{k3}U\left(y\right)=\mathrm{0,}k=\mathrm{1,2,}$ (20)

where $\Delta =\frac{{\partial }^2}{\partial y_1^2}+\frac{{\partial }^2}{\partial y_2^2}$ is the two-dimensional Laplace operator with respect to the local variable $y.$

Substituting Eq. (8) in Eq. (5) yields ${\epsilon }^n{\sigma }_{i3}^{(n)}(x,y)=\mathrm{0,}n=\mathrm{0,1,2...,}i=\mathrm{1,2...,}$and letting $\epsilon \to 0$we obtain the zeroth-order term ${\sigma }_{i3}^{(0)}(x,y)n_i=0$. 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), we obtain ${\sigma }_{i3}^{(0)}(x,y)n_i=\left[c_{44}(y)n_i+{}_{k3}\sigma {}_{i3}(y)n_i\right]\frac{\partial u_{}^{(0)}(x)}{\partial x_k},i=\mathrm{1,2}(nosumoveri).$ Since we have $\frac{\partial u_{}^{(0)}(x)}{\partial x_k}\ne \mathrm{0,}$ we obtain the local boundary equations

${}_{k3}\sigma {}_{i3}(y)n_i=-c_{44}(y)n_k,i,k=\mathrm{1,2.}$ (21)

Substituting Eq. (16) into Eq. (21), we obtain

$\frac{\partial {}_{k3}U{}_{}\left(y\right)}{\partial y_i}{n}_i=-n_k,k=\mathrm{1,2}.$ (22)

For the problem at hand, which is to calculate analytically the effective shear elastic modulus in the longitudinal direction, $c_{44}^{ef},$of the elastic solid, only the Local Problem ${}_{13}L$ over the local domain,$Y,$ illustrated on the right-hand side of Fig. 1, will be specified at the outset and solved. This problem consists of finding a function $U={}_{13}U:Y\to ℝ$ from Eq. (20)-(22), harmonic and of zero average in $Y,$that satisfies the system of equations

$\Delta U=0inY,$ (23)

$U_{\mathrm{,1}}{n}_1+U_{\mathrm{,2}}{n}_2=-n_{1}on\Gamma ,$ (24)

$\langle U\rangle =\mathrm{0,}$ (25)

where $\Delta$ is the two-dimensional Laplacian in $Y,$ $\Gamma$ corresponds to the boundary of the hole in the cell, $F,i=\partial F/\partial y_i$ and $\langle U\rangle \triangleq {\displaystyle \underset{Y}{\int }U\left(y\right)dy/\left|Y\right|,}$with $\left|Y\right|$ being the volume of $Y$. Nevertheless, by solving Local Problem ${}_{23}L$, we verify that the effective shear elastic modulus $c_{55}^{ef}$ is equal to $c_{44}^{ef}$.

The solution of Eq. (23)-Eq. (25) is sought in the class of doubly periodic harmonic functions,$U(z)$, which depends on the complex variable $z=y_1+i{y}_2$, in the form of a series with undetermined real coefficients $a_k$, given by

$U(z)=\mathrm{Re}\{a_1[\zeta (z)-\pi z/\sin \mu ]+{\displaystyle \sum _{k=\mathrm{3,5,...}}^{\infty }{a}_k{\zeta }^{\left(k-1\right)}(z)/(k-1)!}\},$ (26)

where $\mu$ is the angle of the unit cell, which is $\pi /3$ for the hexagonal cell, and $\zeta (z)$ and ${\zeta }^{(k)}(z)$are, respectively, the quasi-periodic Weierstrass zeta function and its kth derivatives of periods ${\omega }_1=1$ and ${\omega }_2=e^{i\mu }$. 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

$c_{44}^{ef}=\langle c_{44}\rangle +c_{44}\langle U{,}_1\rangle ,$ (27)

where $c_{44}$ is the shear modulus of the solid part of the elastic medium and it can be shown (Bravo-Castillero et al., 2009) that

$\langle U{,}_1\rangle =\left(\pi R^2-2\pi a_1\right)/\sin \mu ,$ (28)

where $a_1$ is the solution of the infinite system

$a_l=R^{2l}[\left(1-\pi a_1/\sin \mu \right){\delta }_{1l}+{\displaystyle \sum _{k=\mathrm{1,3,5,...}}^{\infty }{a}_k{\eta }_{kl}}]forl=\mathrm{1,3,...,}$ (29)

where ${\delta }_{1l}$ is the Kronecker delta and

${\eta }_{(kl)}=C_{k+l-1}^l{S}_{k+l},S_{k+l}={\displaystyle \sum _{m,n}*{\beta }_{mn}^{-k-l}}fork+l\ge 2,{\beta }_{mn}^{}=m{\omega }_1+n{\omega }_2,m,n=\mathrm{0,}\pm \mathrm{1,}\pm \mathrm{2,...},m,C_k^l=k!/l!\left(k-l\right)!.$

The star, $*$, in the summation sign indicates that the double sum excludes $m=n=0$. For the hexagonal array of holes considered in this work, the series $S_{k+l}$and ${\eta }_{(kl)}$ are not null if $k+l=6t$ for $t=\mathrm{1,2,3,...}$. The solution $a_1$of the infinite system in Eq. (29) is given by

$a_1=\frac{R^2}{1+A_2-V_p^T{M}_p^{-1}{\tilde{V}}_p},$ (30)

in which $A_2=\pi R^2/\sin \mu$ and we recall from above that $\mu =\pi /3$for a hexagonal cell. Also, for the case of the distribution of hexagonal cells considered in this work,

$V_p^{}=V_p^{}\left({\omega }_s\right),M_p^{}=M_p^{}\left(m_{ts}\right),{\tilde{V}}_p={\tilde{V}}_p\left({\tilde{\omega }}_t\right)$

are vectors and matrices of infinite dimensions and

${\omega }_s=R^{12s}{\eta }_{(16s-1)},m_{ts}={\delta }_{(6t-16s-1)}-R^{12s}{\displaystyle \sum _{i=1}^{\infty }{R}^{12i}{\eta }_{(6t-16i+1)}{\eta }_{(6i+16s-1)},}{\tilde{\omega }}_t={\eta }_{(6t-11)},t,s=\mathrm{1,2,3,...}.f$

Upon substituting Eq. (30) in Eq. (28) and the resulting expression in Eq. (26), we obtain the analytical expression

$c_{44}^{ef}=c_{44}\left(1-2A_2{K}^0\right),$ (31)

where

$K^0=\frac{1}{1+A_2-V_p^T{M}_p^{-1}{\tilde{V}}_p}.$ (32)

4 COMPUTATIONAL METHOD

Let ${\Omega }_H$ be a rectangular region containing a part of the periodic distribution of hexagonal unit cells considered in Section 2 and let ${\Omega }_L$ be a thin layer of thickness $a$ surrounding ${\Omega }_H$, as illustrated by the blue frame in Fig. 2. The layer has the same material properties of the solid part occupying ${\Omega }_H$. The whole region ${\Omega }_H\cup {\Omega }_L$ has dimensions $L\times H$. Let also $M_C$ be the whole region excluding the holes, in analogy to the region $M$ of Section 2. The region $M_C$ 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, ${\partial }_1{M}_C$, zero tractions on the lower, ${\partial }_2{M}_C$, and upper, ${\partial }_4{M}_C$, sides, and displacement $\overline{u}$ on the right-hand side, ${\partial }_3{M}_C$, which has magnitude $\overline{u}$ and is applied normal to the $x_1{x}_2$- plane. Also, $n$ is the outward normal vector to the external boundary.

Figure 2
Cross section of the cylindrical sample together with boundary conditions and the unit cell of hexagonal shape.

To solve numerically the equilibrium problem described at the end of in Section 2 in the domain $M_C$, we observe that this problem is analogous to a linear steady state heat conduction problem, which consists of finding the temperature field $T:M_C\to ℝ$ that satisfies

${\Delta }_{x}T=0in{M}_C,$ (33)

$n\cdot {\nabla }_{x}T=0on{\partial }_H{M}_C\cup {\partial }_2{M}_C\cup {\partial }_4{M}_C,$ (34)

$T=0on{\partial }_1{M}_C,T=\overline{T}on{\partial }_3{M}_C,$ (35)

where ${\Delta }_x$and ${\nabla }_x$ are, respectively, the Laplacian and the gradient operators defined in $M_C$, ${\partial }_H{M}_C$ is the union of all the contours of the holes, and ${\cup }_{i=1}^4{\partial }_i^{}{M}_C$ is the external boundary of $M_C$. We use the finite element commercial package COMSOL 4.4® to obtain approximate solutions of Eq. (33)-Eq. (35). The temperature field $T$is approximated by second-order 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 $\overline{T}$ with the displacement $\overline{u}$ and the conductivity $k$ with the elastic modulus $c_{44}^{}$. The effective elastic modulus $c_{44}^{ef}$ can then be determined from the expression

$c_{44}^{ef}=\frac{{\displaystyle {\int }_{{\partial }_3{M}_C}n\cdot {\nabla }_{x}T}d{x}_2}{\overline{T}/L},$ (36)

where ${\displaystyle {\int }_{{\partial }_3{M}_C}n\cdot {\nabla }_{x}T}d{x}_2$ is the heat flux on ${\partial }_3{M}_C$.

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 $a$ 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 $a$. 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 m2 yields h = 0.0062 m. We then consider both a volume fraction of the holes $A_2=0.4$ 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 $L\times H$ of $M_C$ in Fig. 2 according to the cases (a) $0.0756\times 0.0714m^2$, (b) $0.1623\times 0.1520m^2$, (c) $0.3173\times 0.2947m^2$, (d) $0.6397\times 0.5925m^2$. In Fig. 4 the dimensions of $M_C$ are fixed and given by $L\times H=0.1\times 0.1m^2$, the layer thickness is $a=0.001m$, the volume fraction of the holes is $A_2=0.4$, and we vary the dimension $h$ 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 $M_C$ 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 $M_C$ has no influence on the convergence analysis.

Figure 3
Cross sections of the cylindrical holes according to the FLS approach for$a=0.0031m$, $A_2=0.4$, and $L\times H$with values in (a)$0.0756\times 0.0714m^2$, (b)$0.1623\times 0.1520m^2$, (c)$0.3173\times 0.2947m^2$, (d)$0.6397\times 0.5925m^2$.

Figure 4
Cross sections of the cylindrical holes according to the FDS approach for$L\times H=0.1\times 0.1m^2$, $a=0.001m$, $A_2=0.4$, and $h$with value in (a) h = L/4, (b) h = L/8, (c) h = L/16, (d) h = L/32.

In Section 5 we present results for the effective modulus $c_{44}^{ef}$ obtained from these discretizations. In the first two cases, corresponding to the values $A_2=\mathrm{0.01,}and{A}_2=0.04$, 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 ${\Omega }_L$. Intersections occur for $A_2$ near $0.2$. The case of $A_2=0.4$ is then introduced as an intermediate case to show these intersections. Finally, $A_2=0.9069$ represents the limit case where the walls of the holes are nearly touching each other.

Figure 5
FEM discretizations of the domain$M_C$on the left side and a zoom in of the top left side of these discretizations on the right side. (a) and (b)$A_2=0.01$, (c) and (d) $A_2=0.04$.

Figure 6
FEM discretizations of the domain $M_C$on the left side and a zoom in of the top left side of these discretizations on the right side. (a) and (b)$A_2=0.4$, (c) and (d) $A_2=0.9069$.

5 RESULTS AND DISCUSSION

Here, we present and discuss some analytical and computational results concerning the determination of the effective elastic modulus $c_{44}^{ef}$. 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 $c_{44}^{}=1GPa$ as the elastic modulus of the solid occupying either the region $M$, which has infinite dimensions, or $M_C$, which has finite dimensions.

First, we use the FDS approach, with fixed dimensions of $M_C$ given by $L\times H=0.1\times 0.1m^2$ and fixed layer thickness given by $a=0.0001m$, and consider a sequence of decreasing values of the hexagon dimension $h$ in the set $h=\left\{L/2,L/\mathrm{4,}L/\mathrm{8,}L/\mathrm{16,}L/32\right\}$ Except for h = L/2, the other values of $h$ correspond to the cross sections in Fig. 4. In Fig. 7 we show graphs of the effective modulus $c_{44}^{ef}$ versus the void volume $A_2$ for this sequence. In the figures below, the label $F{i}^{\#}$ stands for a Finite element result obtained with a layer thickness $a$ having $i+1$+1 significant digits and using the FDS approach. Thus, F3# in Fig. 7 means that we have used the FDS approach together with $a=0.0001m$. 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$h\to 0$, even at a high void volume fraction$A_2$. These results show the fact that the numerical result for this constant is influenced by the decreasing values of the hexagon side dimension $h$ 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 $\epsilon$ in Eq. (7), which, here, can be taken as either $h/L$ or $h/H$.

Figure 7
Effective elastic modulus $c_{44}^{ef}$ versus void volume fraction$A_2$using the FDS approach depicted in Fig. 4.

Next, we compare the best numerical curve presented in the Fig. 7, corresponding to $F{3}^{\#}h=L/32$, with the curve obtained via AHM and the best curve by using FLS approach, corresponding to $F4$ obtained by Aguiar et al. (2018). In the figures below, the label $Fi$, 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 $L\times H=0.6397\times 0.5925m^2$ and the thickness of the boundary layer $a=0.0031m$. This figure shows that for a thin boundary layer $a$ 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 $h=L/32$. 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 $h=L/32$ as a representative parameter of number of holes in a fixed domain and below we will vary the boundary layer thickness $a$

Figure 8
Effective elastic modulus $c_{44}^{ef}$ versus void volume fraction$A_2$ using results of Fig. 2 to compare with the corresponding curve$FEMF4$ from Aguiar et al. (2018).

Next, we compare results obtained from both approaches FDS and FLS. In the FDS approach, we still take $L\times H=0.1\times 0.1m^2$, but hold the side of the hexagon fixed at $h=L/32$and consider a sequence of decreasing values of the layer thickness $a.$ in the set $\left\{\mathrm{0.01,}\mathrm{0.001,}\mathrm{0.0001,}0.00001\right\}m$. 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 $c_{44}^{ef}$versus the void volume $A_2$ for both sequences. In Fig. 9(a) we consider the whole interval $\left(\mathrm{0,1}\right)$ and in Fig. 9(b) we show a zoom in of the graphs for the interval$\left(\mathrm{0.75,0.86}\right)$. 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 $Fi{}^{\#},i=\mathrm{1,}\dots \mathrm{,4}$, converge to limit values that are close to the analytical results obtained from AHM as $a$ 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 $\epsilon$ in Eq. (7), which, again, can be taken as either $h/L$ or $h/H$, decreases as both $L$ and $H$increase for a fixed $h$. 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 $A_2=0.9069$ and corresponds to the cases (c) and (d) of Fig. 6, no convergence was obtained.

In Fig. 10 we show magnifications of the upper left corners of the FEM discretizations of $M_C$ having void volume fraction $A_2=0.4$ in the case of the FDS approach. The figure shows how the mesh is modified as the layer thickness $a$ decreases, with values in the set $\left\{\mathrm{0.01,}\mathrm{0.001,}\mathrm{0.0001,}0.00001\right\}m$.

Figure 9
Effective elastic modulus $c_{44}^{ef}$versus void volume fraction$A_2$considering both approaches FDS ($Fi{}^{\#}$) and FLS ($Fi$) and the AHM solution. (a)$A_2\in \left(\mathrm{0,1}\right)$, (b) $A_2\in \left(\mathrm{0.75,0.86}\right)$.

Figure 10
Magnification of the upper left regions of FEM discretizations of $M_C$for$A_2=0.4$considering the FDS approach: (a) $a=0.01m$, (b) $a=0.001m$, (c)$a=0.0001m$, (d) $a=0.00001m$.

In addition to the calculation of the effective modulus $c_{44}^{ef}$ 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 $c_{44}^{ef}$ versus $A_2$ 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 $c_{44}^{ef}$.

Figure 11
Curves of effective elastic modulus $c_{44}^{ef}$ versus void volume fraction $A_2$. Square (SQ), and hexagonal (HEX) lattices together with the FLS (F4) and the FDS (F4#) approaches. (a) $A_2$ in (0,1), (b) $A_2$in (0.005, 0.04).

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 $A_2$, which has the values $0.01$ and $0.04$, the effective shear elastic coefficient obtained by Swan et al. (2003) from the fully drained case, the coefficient $c_{44}^{ef}$ obtained via both AHM and FEM with the FLS approach considered in this work, and percentage errors between their effective coefficients and $c_{44}^{ef}$ obtained via both AHM and FEM. Observe from their Table 2 together with their Table 1 that $A_2=0.01$ is associated with Young’s modulus $E=11GPa$ and Poisson’s ratio $\nu =0.39$and that $A_2=0.04$ is associated with $E=12GPa$ and $\nu =0.38$. 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$A_2=0.04$ and a larger discrepancy in the case$A_2=0.01$, which is being investigated.

Table 1
calculated effective constant and relative error between results obtained via AHM and FEM (FLS approach) and results presented by Swan et al. (2003) for the fully drained case.

$A_2$	(GPa) Swan et al. (2003) Fully drained	(GPa) AHM	(GPa) FEM	Relative Error (%) AHM versus Swan et al. (2003)	Relative Error (%) FEM versus Swan et al. (2003)
0.01	3.5587	3.8784	3.8776	8.9836	8.9613
0.04	4.0161	4.0134	4.0098	0.0672	0.1572

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 $c_{44}^{ef}$ 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.

6 CONCLUSION

In this paper, we have shown some results concerning the analytical and numerical evaluations of the effective modulus $c_{44}^{ef}$ of linear elastic solids having a periodic distribution of voids arranged in hexagonal cells. When the side of the hexagon with length $h$decreases and the boundary layer thickness $a$ is fixed, or, when $a$ decreases and $h$is fixed, we observe convergence of $c_{44}^{ef}$ 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 $c_{44}^{ef}$ that are very close to each other up to $A_2=0.5$ for both $h/L$ and $h/H$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.

Acknowledgments

Financial support is gratefully acknowledged from National Council for Scientific and Technological Development (CNPq), grants # 420099/2018-2, and São Paulo Research Foundation (FAPESP), grant # 2014/21836-2. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001. We would also like to express our gratitude to Prof. Dr. Daniel Varela Magalhães for kindly providing the use of license of COMSOL 4.4® (COMSOL Multiphysics (2013).

References

• Aguiar, A. R., Bravo Castillero, J., Rodríguez Ramos, R., And Da Silva, U. P. (2013). Effective electromechanical properties of 622 piezoelectric medium with unidirectional cylindrical holes. ASME. J. Appl. Mech 80(5): 050905.
• Aguiar, A. R., Prado, E. B. T., Da Silva, U. P. (2018). Effective moduli of 3-1 longitudinally porous solids with regular hexagonal array. In: 6º Encontro Nacional de Engenharia Biomecânica - 6º ENEBI, 2018, Águas de Lindóia. Anais do 6º Encontro Nacional de Engenharia Biomecânica. Rio de Janeiro: ABCM, 2018, 1.
• Bravo-Castillero, J., Rodriguez-Ramos, R., Guinovart-Diaz, R., Sabina, F. J., Aguiar, A. R., Da Silva, U. P., Gomez-Munoz, J. L. (2009). Analytical formulae for electromechanical effective properties of 3-1 longitudinally porous piezoelectric materials. Acta Materialia 57: 796-803.
• COMSOL Multiphysics (2013). Reference Manual, version 4.4. COMSOL, Inc., www.comsol.com
  » www.comsol.com
• Ćwikła, G., Grabowik, C, Kalinowski, K., Paprocka, I., Ociepka, P. (2017). The influence of printing parameters on selected mechanical properties of FDM/FFF 3D-printed parts. IOP Conference Series: Materials Science and Engineering 227: 012033
• Ghosh, S., Lee, K., and Raghavan, P. (2001). A multi-level computational model for multi-scale damage analysis in composite and porous materials. International Journal of Solids and Structures 38(14): 2335-2385.
• Ghosh, S., Lee, K., Moorthy, S. (1996). Two scale analysis of heterogeneous elastic-plastic materials with asymptotic homogenization and Voronoi cell finite element model. Computer Methods in Applied Mechanics and Engineering 132: 63-116.
• Grimal, Q., Raum, K., Gerisch, A., Laugier, P. (2008). Derivation of the mesoscopic elasticity tensor of cortical bone from quantitative impedance images at the micron scale. Computer Methods in Biomechanics and Biomedical Engineering 11 (2): 147-157.
• Hollister, S.J., Kikuchi, N. (1992). A comparison of homogenization and standard mechanics analyses for periodic porous composites. Computational Mechanics 10: 73-95.
• Lanzotti, A., Grasso, M., Staiano, G., Martorelli, M. (2015). The impact of process parameters on mechanical properties of parts fabricated in PLA with an open-source 3-D printer. Rapid Prototyping Journal 21 (5): 604-617.
• Lo Cascio, M., Milazzo, A., Benedetti, I. (2020). Virtual element method for computational homogenization of composite and heterogeneous materials. Composite Structures 232: art. no. 111523.
• López-Realpozo, C., Rodríguez-Ramos, R., Guinovart-Diaz, R., Bravo-Castillero, J., Sabina, F.J. (2011). Transport properties in fibrous elastic rhombic composite with imperfect contact condition. International Journal of Mechanical Sciences 53: 98-107.
• López-Realpozo, C., Rodríguez-Ramos, R., Guinovart-Diaz, R., Bravo-Castillero, Otero, J.A., Sabina, F.J., Lebon, F., Dumont, S., Sevostianov, I. (2014). Effective elastic shear stiffness of a periodic fibrous composite with non-uniform imperfect contact between the matrix and the fibers. International Journal of Solids and Structures 51: 1253-1262.
• Marino, M., Hudobivnik, B., Wriggers, P. (2019). Computational homogenization of polycrystalline materials with the Virtual Element Method. Computer Methods in Applied Mechanics and Engineering 355: 349-372.
• Otero, J. A., Rodríguez-Ramos, Bravo-Castillero, J., Guinovart-Díaz, R., Sabina, J., Monsivais, G. (2013). Semi-analytical method for computing effective properties in elastic composite under imperfect contact. International Journal of Solids and Structures 50: 609-622.
• Otero, J. A., Rodríguez-Ramos, R., Monsivais, G. (2016).Computation of effective properties in elastic composites with different inclusion shapes and under imperfect contact. Mathematical Methods in the Applied Sciences 40 (9): 3290-3310.
• Parnell, W., Grimal, Q. (2009). The influence of mesoscopic porosity on cortical bone anisotropy. Investigations via asymptotic homogenization. J R Soc Interface 6:97-109.
• Parnell, W.J., Vu, M.B., Grimal, Q., Naili, S. (2012). Analytical methods to determine the effective mesoscopic and macroscopic elastic properties of cortical bone. Biomech Model Mechanobiol 11:883-901.
• Pingaro, M., Reccia, E., Trovalusci, P. (2019). Homogenization of Random Porous Materials With Low-Order Virtual Elements ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems, Part B: Mechanical Engineering, 5 (3), art. no. 030905.
• Provaggi, E., Capelli, C., Rahmani, B., Burriesci, G, Kalaskar, D. M. (2019). 3D printing assisted finite element analysis for optimising the manufacturing parameters of a lumbar fusion cage, Materials & Design 163: 107540.
• Ramakrishnan, N., and Arunachalam, V. S. (1990). Effective elastic moduli of porous solids. Journal of Materials Science 25(9): 3930-3937.
• Ramírez-Torres, A., Penta, R, Rodríguez-Ramos, R., Merodio, J., Sabina, F. J., Bravo-Castillero, J, Guinovart-Díaz, R., Preziosi, L., Grillo, A. (2018). Three scales asymptotic homogenization and its application to layered hierarchical hard tissues. International Journal of Solids and Structures 130-131: 190-198.
• Said, B. M., Salah, M., Kanit, T., and Kamel, F. (2016). On the homogenization of 2d porous material with determination of rve. International Journal of Mechanical and Mechatronics Engineering 16(1): 81-86.
• Sevostianov, I., Kachanov, M. (2000). Impact of the porous microstructure on the overall elastic properties of the osteonal cortical bone. ASME. J Biomech 33(7):881-888.
• Swan, C.C., Lakes, R.S., Brand, R.A., Stewart, K.J. (2003). Micromechanically based poroelastic modeling of fluid flow in haversian bone. ASME. J Biomech Eng 125(1):25-37.
• Wohlers, T. (Ed.) (2012). Wohlers Report 2012: additive manufacturing and 3D printing state of the industry. Wohlers Associates Inc.
• Yap, Y. L., Toh, W., Koneru, R., Lin, K., Yeoh, K. M., Lim, C. M., Lee, J. S., Plemping, N. A., Lin, R., Ng, T. Y., Chan, K. I., Guang, H., Chan, W. Y. B., Teong, S. S., Zheng, G. (2019). A non-destructive experimental-cum-numerical methodology for the characterization of 3D-printed materials-polycarbonate-acrylonitrile butadiene styrene (PC-ABS). Mechanics of Materials 132: 121-133.

Edited by

Guest Editors:

Volnei Tita and Nicholas Fantuzzi.

Publication Dates

• Publication in this collection
  13 Nov 2020
• Date of issue
  2020

History

• Received
  31 Jan 2020
• Reviewed
  07 Aug 2020
• Accepted
  28 Aug 2020
• Published
  08 Sept 2020