Scale-size effects analysis of optimal periodic material microstructures designed by the inverse homogenization method

Computers and Structures xxx (2015) xxx–xxx

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Scale-size effects analysis of optimal periodic material microstructures designed by the inverse homogenization method P.G. Coelho a,⇑, L.D. Amiano b, J.M. Guedes c, H.C. Rodrigues c a

UNIDEMI, DEMI, Faculty of Sciences and Technology, Universidade Nova de Lisboa, Portugal DEMI, Faculty of Sciences and Technology, Universidade Nova de Lisboa, Portugal c IDMEC, Universidade de Lisboa, Portugal b

a r t i c l e

i n f o

Article history: Accepted 7 October 2015 Available online xxxx Keywords: Homogenization Microstructures Optimization Topology Cellular Scale

a b s t r a c t Periodic homogenization models are often used to compute the elastic properties of periodic composite materials based on the shape/periodicity of a given material unit-cell. Conversely, in material design, the unit-cell is not known a priori, and the goal is to design it to attain specific properties values – inverse homogenization problem. This problem is often solved by formulating it as an optimization problem. In this approach, the unit-cell is assumed infinitely small and with periodic boundary conditions. However, in practice, the composite material comprises a finite number of measurable unit-cells and the stress/strain fields are not periodic near the structure boundary. It is thus critical to investigate whether the obtained topologies are affected when applied in the context of real composites. This is done here by scaling the unit-cell an increasing number of times and accessing the apparent properties of the resulting composite by means of standard numerical experiments. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction The optimal design of materials with periodic microstructure to achieve specific properties has been a major area of interest in materials research [1]. The so-called ‘‘unit-cell” is the microstructure representative of the smallest periodic heterogeneity of the material domain. The material is then defined assembling unitcells, i.e. repeating the unit-cell in all directions of the space. Periodic cellular composites are so attractive because they offer great flexibility in tailoring different physical properties by controlling the microstructural architecture in unit-cells and they can be realized by additive manufacturing technologies [2]. This has opened the way to explore and develop novel multifunctional composites even to the point of exhibiting exceptional properties usually not encountered in nature – metamaterials [3]. To take full advantage of the possibilities provided by these materials systematic design methods are required. The inverse homogenization method is often used for this purpose [4]. The goal is to design the unit-cell targeting prescribed effective (homogenized) properties or extremizing a function of those properties ⇑ Corresponding author. E-mail addresses: [email protected] (P.G. Coelho), [email protected] (L.D. Amiano), [email protected] (J.M. Guedes), [email protected] (H.C. Rodrigues).

(e.g. strain energy). This can be solved by designing the unit-cell through topology optimization [5]. Material type or density in each element of the finite element mesh discretizing the unit-cell serves as the design variable and hence, the design freedom is huge allowing for exceptional improvements on material efficiency. A thorough review of periodic microstructure design using inverse homogenization supported on topology optimization was carried out by Cadman et al. in 2013 [1]. Further advances in the design of periodic microstructures of multi-functional materials for a range of physical properties have been published in the meantime which proves how active this research area has been. For instance, quite recent research papers present optimal material designs for poroelastic actuators [6], elastic stiffness [7–11], Poisson’s ratio [2,9], natural frequencies [12], thermal conductivity [10], fluid permeability [11], negative or zero compressibility [13], viscoelastic behavior [14–16], bone scaffolds [17], piezocomposites [18,19]. The inverse homogenization method assumes that the scale of a unit-cell is indeterminate (infinitely small) as well as periodic Boundary Conditions (BC’s) [20–23]. This makes uncertain whether the obtained topology can be translated into real composites of macroscale. In practice, one has a finite number of measurable unit-cells assembled together to define the composite sample. Moreover, the stress or strain fields are in general arbitrary (not periodic) near the boundary of the composite.

http://dx.doi.org/10.1016/j.compstruc.2015.10.001 0045-7949/Ó 2015 Elsevier Ltd. All rights reserved.

Please cite this article in press as: Coelho PG et al. Scale-size effects analysis of optimal periodic material microstructures designed by the inverse homogenization method. Comput Struct (2015), http://dx.doi.org/10.1016/j.compstruc.2015.10.001

2

P.G. Coelho et al. / Computers and Structures xxx (2015) xxx–xxx

Therefore it is critical to investigate how the apparent properties of the composite converge as one scales the unit-cell an increasing number of times. Numerical standard experiments based on BC’s of uniform stress and strain can be applied to access the apparent properties of the resulting composite for a specific scale factor [24]. This factor is defined here as the ratio between the linear size of the composite and the linear size of the unitcell. As this factor increases one may evaluate convergence of the apparent properties to the theoretical values from homogenization. This scale-size effects analysis has been carried out by other authors although restricted to material symmetries, mostly to two-dimensional microstructures and investigating mainly the convergence of the mean compliance, few elastic moduli and fundamental frequency [25–32]. The aim of the present work is to extend this to three-dimensional anisotropic material microstructures. This paper is based upon Coelho et al. [33], but the current paper adds the full apparent elasticity tensor estimate based on additional numerical experiments, measures the degree of microstructural anisotropy, investigates convergence on more elastic tensor coefficients as well as on strain/stress energy density and increases also the scale factor from 5 to 6. Ideally a convergence analysis of optimal periodic microstructures would be carried out, from the finite periodicity to the theoretically infinite periodicity. This extreme case is of course practically infeasible. Due to the limitation of the computational resources available, only a few low scale factors have been tried for three-dimensional analyses. Despite this analysis limitation, the outcome of this study indicates that it is sufficient to have a low scale factor to replace the non-homogeneous composite by the equivalent homogeneous material with the moduli computed by homogenization. This is demonstrated here with material microstructures that extremize elastic stiffness considering different load cases and subjected to either volume fraction or permeability constraints. For the unit-cell topologies selected, one studies ultimately how rapidly the 21 independent elastic coefficients from the apparent fourth-order stiffness tensor converge to the corresponding values of the homogenized stiffness tensor. 2. Materials and methods

solid

Fig. 1. Material model of solid-void composite. Discretization of Y with periodic BC’s. Array of 5  55 unit-cells of global size D (in theory goes to infinity) and one unit-cell of size d.

2.2. Optimization problem In this work one extremizes an energy density based objective function subjected to a minimum prescribed permeability of the material microstructure in all directions of the space. Such a problem is of practical interest, for example, to design scaffolds for skeletal Tissue Engineering. Ensuring isotropic permeability in scaffolds favors the migration of cells, nutrients and waste products needed for native tissue regeneration while ensuring a preferential stiffness direction helps the scaffold fulfilling its rule as a load bearing device. An appropriate trade-off between these two conflicting properties, permeability (favors porosity) and stiffness (favors density), can be achieved through topology optimization [11,34–39]. The local anisotropic material design problem can be then defined by Eq. (1) and a brief description is given below (for full details see e.g. [33,35,40,41]).

1 H ðlÞrmn rkl C 2 mnkl subject to: min

2.1. Material model The material model here is a three-dimensional porous composite material (solid-void phases) generated by the repetition of a unit-cell in all directions of the space. The unit-cell represents the smallest periodic heterogeneity of the composite domain W whose material properties are calculated by homogenization. This theory assumes periodic boundary conditions (BC’s) applied to the unit-cell domain Y and infinite periodicity of the unit-cell (Y-periodic), i.e. its feature size d is much smaller than the cellular material global size D (d/D ! 0), see Fig. 1. In the case of optimal material design, the design domain is the unit-cell domain Y. Ideally, one would determine the pointwise material distribution in this domain that optimizes the design goal. In practice, the true pointwise material distribution cannot be found, but one can discretize the domain Y with a large number of finite elements and let the density l in each element be a variable. Here one uses a regular mesh 20  20  20 of 8-node isoparametric hexahedral finite elements as shown in Fig. 1. Solid and void phases correspond to l equal to 1 and 0, respectively. Topology optimization can then be used to search for an optimal material microstructure layout targeting desired effective properties. This design method is also known as inverse homogenization and it was originally introduced by Sigmund [4]. The effective material properties (homogenized) can be found using a numerical homogenization procedure described in [23].

l

K Hij ðlÞ P k ; 

K Hij ðlÞ ¼ 0;

i ¼ j ¼ 1; . . . ; 3

ð1Þ

i – j and i; j ¼ 1; . . . ; 3

In Eq. (1) l is the local density varying between 0 (void) and 1 (material) which depends on the spatial variable y in the unit-cell design domain Y (see Fig. 1). The problem is formulated for the minimization of the stress energy density or mean compliance. The stress tensor is r and characterizes an averaged macroscopic stress field applied to the composite which is here fixed a priori (e.g. hydrostatic). The homogenized compliance tensor is CH , i.e. the inverse of the homogenized stiffness tensor EH . The stiffness tensor at each point of Y is related to the tensor E0 of the base material properties through the SIMP interpolation scheme [42]. Regarding the optimization problem constraints, one enforces the permeability tensor to be diagonal and each diagonal coefficient has to be  ‘‘equal” or ‘‘greater than” a threshold, k . This way one gets an interconnected pore network on the periodic material satisfying a critical (minimum) permeability in all direction of the space. The permeability measure of porous media is given by the homogenized tensor KH . This tensor is obtained homogenizing a potential flow problem in periodical porous media characterized by the Darcy law [20,22]. Here one considers the interpolation between permeability and local density l given by a power-law

Please cite this article in press as: Coelho PG et al. Scale-size effects analysis of optimal periodic material microstructures designed by the inverse homogenization method. Comput Struct (2015), http://dx.doi.org/10.1016/j.compstruc.2015.10.001

3

P.G. Coelho et al. / Computers and Structures xxx (2015) xxx–xxx

(see e.g. [33–35]). The permeability tensor for the base material, K0 , is here considered to be unitary, diagonal and isotropic. This way the interpolation scheme means that void and solid have high (100%) and low (0%) permeability, respectively. Consequently, the homogenized permeability tensor KH characterizing the periodic material medium becomes normalized, i.e. its coefficients take val ues between 0 and 1 (0% and 100%), and thus the threshold k in Eq. (1) can also be fixed between these bounds. One uses the adjoint sensitivity analysis [5,33] to obtain the gradients and the Method of Moving Asymptotes [43] for the design updates. Finally, one remarks that imposing different  thresholds k in Eq. (1) one attains different achievable ranges on solid volume fraction, q, defined by Eq. (2). This way one may approach the actual bound that relates permeability and solid volume fraction [44,45].

q¼

Z

Y

lðyÞdY:

ð2Þ

2.3. Measuring anisotropy The compliance tensor CH can be rotated using a rotation tensor R in terms of the Euler angles, according to Eq. (3). At each rotation the component 1=C 0H 1111 has the meaning of a directional stiffness and it is used here for plotting as a measure of the degree of material anisotropy. In the particular case of isotropy, that quantity is the Young’s Modulus E and the corresponding plot is a sphere (see also [46]). H C 0H ijkl ¼ Rim Rjn Rkp Rlq C mnpq

ð3Þ

As an addition to the graphical representation of material anisotropy, one also measures anisotropy in terms of a scalar A defined in Eq. (4), see also [47]. This is a measure interpreted as the ‘‘distance” between an effective elasticity tensor Eijkl and the ‘‘nearest” isotropic elasticity tensor E;iso ijkl . Its value ranges from 0 (isotropy) to 2 (e.g. unidirectional stiffness as shown later in Section 3).

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   u  u E  E;iso E  E;iso u ijkl ijkl ijkl ijkl A¼t ;iso E;iso ijkl Eijkl

increasing scale factor n. This factor is defined here as the ratio between the linear size of the composite D and the linear size of the unit-cell d; n ¼ D/d, see Fig. 1 and [31]. Assuming that the composite domain W is always here an unitary cube, jWj = 1 or D = 1, n becomes equal to the number of times that the unit-cell is repeated along each Cartesian axis and thus one combines n3 of such cells inside W. Recall that the homogenized elastic properties EH ijkl in Eq. (1) are calculated when n ! 1. Consequently, in the homogenization result there is no length for the unit-cell (in the limit d = 0) and also there is no sensitivity to the variation of the scale factor n. In engineering practice, size d is measurable and the mechanical response of the composite under arbitrary BC’s is n-dependent. That is why, it is critical to investigate how accurate homogenization predictions are compared to the actual properties of a composite with the unit-cell scaled n times. In this work the scale factor n varies from 1 to 6 as illustrated in Fig. 2. Realizing the successive numerical models for the composites in Fig. 2 (from the left to the right) one quickly notices that a prohibitively large computational problem in the span of a few low scale factors can be created. For example, in case of n = 5, the 8000 elements in Fig. 1 multiply by 53 cells and thus the numerical model for the resulting composite has already 1 million elements and more than 3 million degrees of freedom. The maximum scale factor achievable by the computational resources available for the present study was n = 6. In order to estimate apparent elastic properties, this work follows a standard numerical testing procedure based on constant strain or stress as detailed below in Section 2.4.1 (see also [24]). The use of alternative BC’s is also justified and considered later in Section 2.4.2. 2.4.1. Constant strain or stress For instance, let’s consider the Dirichlet-type displacement boundary conditions,

uðxÞj@ W ¼ HðiÞ xj@ W

ð4Þ

rðxÞnj@W ¼ HðiÞ nj@W i ¼ 1 to 6 where

tive tensor Eijkl , i.e. quantities k and l in Eq. (5) that are referred as the bulk and shear modulus in the isotropic case, respectively. 



1  E ; 9 iijj

l ¼

ð6Þ

as well as the Neumman-type traction boundary conditions:

The isotropic tensor E;iso ijkl is defined from the invariants of the effec-

k ¼

i ¼ 1 to 6

 1 1   Eijij þ Eijji  E 20 30 iijj

ð5Þ

2.4. Numerical experiments Numerical testing procedures applying specific BC’s are used here to measure mechanical responses of the composite in order to investigate how the apparent properties converge with an

2

b 6 1 H ¼ 40 0 2 0 6 H4 ¼ 4 b

0 0

3

7 0 0 5; 0 0

3 b 0 7 0 0 5; 0 0 0

2

0 6 2 H ¼ 40 0 2 0 6 H5 ¼ 4 0

ð7Þ

0 0 b

3

7 0 5;

0 0

3 0 0 7 0 b 5; 0 b 0

2

0 6 3 H ¼ 40 0 2 0 6 H6 ¼ 4 0 b

0 0

3

7 0 0 5; 0 b

3 0 b 7 0 05 0 0

ð8Þ

and @ W is the boundary of the composite sample W, u is the displacement vector, x is the spatial position vector, r is the Cauchy stress tensor, n is the outward normal unit vector and b is a constant. The upper script index i on tensor H means the application

Fig. 2. Unit-cube containing arrays n  n  n of unit-cells (n varies from 1 to 6).

Please cite this article in press as: Coelho PG et al. Scale-size effects analysis of optimal periodic material microstructures designed by the inverse homogenization method. Comput Struct (2015), http://dx.doi.org/10.1016/j.compstruc.2015.10.001

4

P.G. Coelho et al. / Computers and Structures xxx (2015) xxx–xxx

of six fundamentals tests (three normal and three shear mechanical tests shown in Eq. 8). In the case of Dirichlet-type conditions, where the applied displacement field is x linearly dependent, the sample W is tested at a uniform strain b. In case of Neumann’s, the test is carried out on W at a uniform stress b. Fig. 3 shows schematically the application of both types of BC’s in a two-dimensional composite domain and one unit-cell (n = 1, but it could be more as indicated in Fig. 2) for the sake of simplicity. Fig. 3a and b illustrates, respectively, the normal (for y-direction) and shear (for in-plane xy) test procedures through the application of a linear displacement field according to Eq. (6). Fig. 3c and d shows, for the same normal and shear directions, the application of the stress-based BC’s according to Eq. (7). The specific nature of the Dirichlet and Neumann-type conditions allow direct estimation of the stiffness and compliance tensors, herein after referred to as Dirichlet ED and Neumann CN tensors, respectively. In both cases, respective average stress hri and strain hei fields are computed through Eq. (9) over all of the N elements in the W mesh, We with e = 1. . .N (both solid and void phases are modelled), and for each one of the six fundamental tests (Eq. 8).

R

hiW ¼

dw jWj

W

ð9Þ

This Eq. (9) is a global average in the sense that is extended to the entire composite volume jWj. Replacing jWj by the volume of the unit-cell jYj, located at the center of W, one obtains a local average based on stress and strain distributions less affected by the arbitrary (non-periodic) BC’s. In this work both averages are taken into account but the data in charts of Section 3 correspond to predictions based on global averages. Once these averages are computed, the Dirichlet and Neumann tensor coefficients can be easily obtained from the material constitutive law, r ¼ Ee, i.e.

hri ¼ ED hei

ð10Þ

hei ¼ C N hri

ð11Þ 

Defining EN ¼ CN

1

as the Neumann stiffness tensor, one replaces

Eq. (11) by

hri ¼ EN hei

ð12Þ

To sum up, on the one hand, one computes Dirichet coefficients EDijkl through Eqs. (6) and (10). On the other hand, one may compute Neumann coefficients ENijkl through Eqs. (7) and (12). In both cases,

one computes 36 constitutive constants Eijkl in the following relation between averages,

(a)

(b)

8 9 2  hr11 iW > E1111 > > > > > 6  > > > hr22 iW > E2211 > > 6 > > > > 6 < hr33 iW = 6 E3311 ¼6 6  > > hs12 iW > > 6 E1211 > > > > 6  > > > > hs23 iW > > 4 E2311 > > : ; hs13 iW E1311

9 38 E1122 E1133 E1112 E1123 E1113 > he11 iW > > > > > he22 i > > E2222 E2233 E2212 E2223 E2213 7 > > 7> W> > > > 7>      < E3322 E3333 E3312 E3323 E3313 7 he33 iW = 7      7 E1222 E1233 E1212 E1223 E1213 7> > hc12 iW > > > > > 7> > > E2322 E2333 E2312 E2323 E2313 5> > hc23 iW > > > > : ; hc13 iW E1322 E1333 E1312 E1323 E1313 ð13Þ

and then these estimates can be directly compared with the homogenized stiffness coefficients EHijkl for each scale factor n. In Eq. (13) only 21 out of 36 coefficients (truly out of 81 coefficients because it’s a fourth-order tensor) are actually independent due to the following symmetries, Eijkm ¼ Ejikm ¼ Eijmk ¼ Ekmij . Nevertheless, the six fundamental tests in Eq. (8) can provide the full characterization of the elastic tensor, as seen in Eq. (13), and one takes advantage of that to check the expected tensor symmetry and thus correctness of the applied BC’s. 2.4.2. Alternative boundary conditions The motivation for considering additional BC’s here has to do with the fact that the Neumann-type conditions underestimate too much the apparent properties of porous materials. In Fig. 3c can be seen the application of the constant pressure b directly on the top of the finite elements with a larger compliance (void) phase. This leads to very high deformations in this phase near the boundary. Thus the elastic coefficient estimative based on Eq. (12) becomes too small since the average strains hei are much higher and the average stress hri remains constant. This unwanted result is a direct consequence of having the void phase here. This wouldn’t be an issue in case of composites mixing two engineering materials (e.g. E-Glass and epoxy resin) where the stiffness ratio could be around 10, see e.g. [29,31]. Keeping the focus of the present work mainly on the analysis of solid-void composites, one overcomes the excessive compliance of the void phase by considering alternative (or mixed) BC’s as summarized in Fig. 4. On the one hand, one pursues a battery of fundamental tests applying Eq. (7) again but the constant pressure is now applied only on the solid phase as illustrated in Fig. 4a and b. Limiting here the analysis to orthotropic material designs, the surface area of the solid phase comes equal between opposite facets but it can differ between adjacent facets. This justifies different applied stresses, s and s0 as illustrated in Fig. 4b, in order to produce equal resultant forces to satisfy Statics. On the other hand, as an alternative to the Neumann-type condition, one applies a uniform pressure upon a rigid plate (here 104 times stiffer than the solid phase) that in turn transfers the load to the top surface of the ‘‘solid” and ‘‘void” elements of the sample

(c)

(d)

Fig. 3. Numerical tests (2D representation of w with n = 1). Dirichlet-type conditions: (a) normal test; (b) shear test. Neumann-type conditions; (c) normal test; (d) shear test.

Please cite this article in press as: Coelho PG et al. Scale-size effects analysis of optimal periodic material microstructures designed by the inverse homogenization method. Comput Struct (2015), http://dx.doi.org/10.1016/j.compstruc.2015.10.001

5

P.G. Coelho et al. / Computers and Structures xxx (2015) xxx–xxx

(a)

(b)

(d)

(c)

Fig. 4. BC’s alternative to Neumann’s to test porous materials (2D representation of w with n = 1): (a) normal traction on solid; (b) shear stress on solid; (c) normal traction upon rigid plate; (d) shear based on displacement BC’s.

preventing this way excessive deformation. One realizes that the Neumann’s assumption of constant traction on the boundary @ W is now no longer satisfied. The plate is modelled with shell elements and it is node-coupled (through the normal direction only) with the mesh of W at the surface. The opposite surface has the displacement constrained in the loading direction, see Fig. 4c. This might be viewed as the equivalent to the test performed on the MTS machine for Modulus assessment (neglecting friction effects). Using rigid plates attached to several facets for shear tests is cumbersome. Hence one decides here to apply displacement-based BC’s. The displacement field given by Eq. (6), i = 4 to 6, is applied taking into account the following simplifications. Displacements are not applied to all @ W but only in two sets of opposite facets and only the displacements that are x-dependent are applied (compare Fig. 3b with Fig. 4d). This is less restrictive than Dirichlet BC’s. Finally, experiments as just described in this subsection are hardly applied to composite designs with arbitrary anisotropy because the satisfaction of static equilibrium equations becomes a critical issue as well as getting symmetry on the effective tensor E in Eq. (13). Therefore, this justifies here the application of these alternative BC’s only to orthotropic designs. 2.5. Energy convergence Several studies on scale-size effects show the mean compliance (energy) convergence [25,26,29–31]. For example, according to the theoretical discussion presented in [31], one expects the chain of inequalities in Eq. (14) to hold. ð1Þ

ð2Þ

ð3Þ

ð1Þ

EN 6 EN 6 EN 6   6 EN

6 E 6 ED

ð1Þ

ð3Þ

ð2Þ

ð1Þ

6   6 ED 6 ED 6 ED

strain energy densities eED e to converge from above while eEN e con-

verge from below to the homogenization result, eEH e, as n increases. Alternatively to Eq. (14), since C ¼ E1 , it holds comprehensively the following chain of inequalities, ð1Þ

ð2Þ

ð3Þ

CD 6 CD 6 CD 6   6 CD

ð1Þ

ð1Þ

6 C 6 CN

ð3Þ

ð2Þ

6   6 CN 6 CN 6 CN

ð1Þ

ð17Þ to be interpreted again in the light of Eq. (15) but now the arbitrary tensor e is made equal to the given macroscopic stress tensor r in Eq (1). Consequently, one expects now stress energy densities

rCN r to converge from above while rCD r converge from below to the homogenization result, rCH r ¼ eEH e, as n increases. The inequalities in Eqs. (14) and (17) are checked in this work. 3. Results Firstly, one shows how the different optimal anisotropic designs obtained approach the bound that relates isotropic permeability with volume fraction. The problem in Eq. (1) is solved first using the material model in Fig. 1 and then the volume fraction of the solid phase, q, is calculated through Eq. (2). As a base material for the unit-cell, one selects polycaprolactone (PCL), often used in biomedical applications, with Young Modulus E = 290 MPa and Poisson coefficient m ¼ 0:3. To the void phase (see Fig. 1) is assigned a comparatively much lower stiffness such that Esolid =Ev oid ¼ 1012 . Imposing different permeability thresholds k in Eq. (1) one attains different achievable ranges on q. Fig. 5 shows this result considering unit-cell designs optimized for different 

ð14Þ Here the superscript denotes the Neumann or Dirichlet tensors of the composite with the scale factor nðn ¼ 1; 2; 3; . . . ; 1Þ. Eq. (14) is not a general result that compares individual tensor coefficients as the scale factor n changes. To interpret correctly these inequalities one remarks that for two general elastic tensors A and B,

B  A P 0 is equivalent to e : ðB  AÞ : e P 0;

8e;

ð15Þ

where e is an arbitrary symmetric second-order tensor. For instance, e can be made equal to the macroscopic strain tensor e evaluated from a given macroscopic stress state r through

e ¼ CH ðl Þr

ð16Þ

where CH ðl Þ is the homogenized compliance tensor evaluated for the optimal design l solution of the problem in Eq. (1). This way Eq. (14) compares strain energy densities. Therefore, one expects

Fig. 5. Permeability threshold against solid volume fraction considering optimization for different load cases. The microstructures shown correspond to the triaxial load case. Results are compared to the Hashin & Shtrikman bound.

Please cite this article in press as: Coelho PG et al. Scale-size effects analysis of optimal periodic material microstructures designed by the inverse homogenization method. Comput Struct (2015), http://dx.doi.org/10.1016/j.compstruc.2015.10.001

6

r and corresponding optimal designs (solid part only) for volume fraction and permeability constraints. Graphical representation of anisotropy and corresponding A.

Stress state (r11 , r22 , r33 , r12 , r23 , r13 Þ

Triaxial

Biaxial

Shear-2

Shear-3

σ

2σ

x2 x3

Uniaxial

1.5σ x1

σ

σ

σ 2σ

3σ

σ

2σ

1.75σ

(r; 2r; 3r; 0; 0; 0)

(0; r; 2r; 0; 0; 0)

(0; 0; r; 0; 0; 0)

(0; 0; 0; r; 1.75r; 0)

(0; 0; 0; r; 2r; 1.5r)

A ¼ 0:3997

A ¼ 0:9359

A ¼ 0:3649

A ¼ 0:7291

A ¼ 0:4200

Volume constraint only V  ¼ 50%

x 10

1

1

0

0

0

-1

y

1

-1 -5

x 10

7

0

5 x

1

0 z

-1 x 10

8

8

x 10

-1 8

0 x

1

1

0 z

-1 x 10

8

1

1 0

x 10

8

0

-1

-1

x 10

8

y

8

y

x 10

y

y

x 10

8

-1 x 10

8

0 x

1

1

0 z

-1 8 x 10

x 10

8

1 x

0 -1

-1 -1 -1

0 z

1

x 10

8

x 10

8

0 x

1 1

0 z

-1 8 x 10



Permeability constraint only k ¼ 50%

A ¼ 1:0806

7

x 10

0

7

0 x

5

5

0 z

-5 7 x 10

-5 x 10

x 10

0 -5

-5

-5 -5

A ¼ 0:3974

7

5 y

0

x 10

x 10

5 y

y

5

A ¼ 2:0000

7

7

0 x

5

5

0 z

-5 7 x 10

A ¼ 0:3790

7

4 2 0 -2 -4 -4 -2 7 0 2 x 10 4 x

x 10

-5 x 10

7

0 x

5

5

0 z

-5

x 10

7

7

4 2 0 -2 -4 -4 -2 7 02 x 10 4 x

y

x 10

y

A ¼ 0:4658

-4 0 -2 7 x 10 4 2 z

-4 0 -2 7 x 10 4 2 z

P.G. Coelho et al. / Computers and Structures xxx (2015) xxx–xxx

Please cite this article in press as: Coelho PG et al. Scale-size effects analysis of optimal periodic material microstructures designed by the inverse homogenization method. Comput Struct (2015), http://dx.doi.org/10.1016/j.compstruc.2015.10.001

Table 1 Macroscopic stress fields

7

P.G. Coelho et al. / Computers and Structures xxx (2015) xxx–xxx Table 2 Homogenized stiffness tensors EH assuming Esolid =Ev oid ¼ 1012 for six microstructural designs selected from Table 1 (units in [MPa]).

-15

E1313 E1313

E3333

E2233

-50

1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6

E2222

200 175 150 125 100 75 50 25 0 -25

E1133

E1313

E2323

E1212

E3333

E2233

E2222

1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6

E2323

(b)

E2323

(a)

E1212

Elastic coeficients

Deviation [%]

Elastic coeficients

E1212

E1313

E2323

E1212

E3333

E2233

E2222

E1133

E3333

-25

-35

E1133

3 0:49 0:13 7 7 0:19 7 7 0:27 7 7 0:27 5 6:64

1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6

-20

-30

E1122

0:14 0:48 0:34 0:21 7:22

-10

0

E1122

1:22 1:08 0:14 6:09

-5

50 E1111

3 0:00 0:00 7 7 0:00 7 7 0:00 7 7 0:00 5 0:00

0

150 100

250 225 200 175 150 125 100 75 50 25 0 -25 -50 -75

0:00 0:00 0:00 0:00 0:00



E2233

250 200

0:00 0:00 0:00 0:00

Shear-3 ðk ¼ 50%Þ 2 43:24 5:22 5:84 6 46:55 6:57 6 6 48:09 6 6 6 4 SYM

3 14:60 4:73 7 7 11:68 7 7 1:47 7 7 13:90 5 33:54

E2222

300

4:02 3:02 1:92 8:00 31:71

E1133

400 350

20:78 18:94 10:11 36:57



Uniaxial ðk ¼ 50%Þ 2 0:00 0:00 0:00 6 0:00 0:00 6 6 81:36 6 6 6 4 SYM

3 0:00 0:00 7 7 0:00 7 7 0:00 7 7 0:00 5 0:00

5

1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6

E1111

Deviation [%]

500 450

Deviation [%]

Shear-3 ðV  ¼ 50%Þ 2 68:95 38:64 35:11 6 99:81 35:44 6 6 122:09 6 6 6 4 SYM

3 25:29 17:57 7 7 33:27 7 7 0:00 7 7 0:00 5 28:45

0:00 0:00 0:00 0:00 5:41

E1122

0:00 0:00 0:00 18:25 39:75

0:00 0:00 0:00 0:00

E1122

0:00 0:00 0:00 19:71



Biaxial ðk ¼ 50%Þ 2 0:00 0:00 0:00 6 49:32 7:48 6 6 67:65 6 6 6 4 SYM

3 0:00 0:00 7 7 0:00 7 7 0:00 7 7 0:00 5 4:72

E1111

0:00 0:00 0:00 0:00 5:89

Deviation [%]

Shear-2 ðV  ¼ 50%Þ 2 24:60 15:17 25:97 6 157:42 35:91 6 6 93:73 6 6 6 4 SYM

0:00 0:00 0:00 5:35

E1111



Triaxial ðk ¼ 50%Þ 2 37:00 5:25 5:24 6 46:80 7:31 6 6 57:00 6 6 6 4 SYM

Elastic coeficients

Elastic coeficients

(c)

(d)

Fig. 6. Scale-size effects analysis for the triaxial case (k ¼ 50%). (a) Dirichlet-type BC’s with Esolid =Ev oid ¼ 1012 ; (b) Neumann-type BC’s with Esolid =Ev oid ¼ 10 ; (c) Constant stress on solid with Esolid =Ev oid ¼ 1012 ; (d) Mixed BC’s with Esolid =Ev oid ¼ 1012 . 

given macroscopic stress fields r pictured in Table 1. The curves shown for the different load cases almost overlap into one located slightly below the theoretical bound from Hashin & Shtrikman [44]. In all optimizations the initial design is the same, i.e. l ¼ 0:5 on vertexes and edges of Y and l ¼ 0:02 in the remaining volume.  Table 1 summarizes the different designs obtained for k ¼ 50% which are basically open wall cells achieving approximately the same volume fraction (around 30%). Replacing the permeability constraint in Eq. (1) by the volume fraction constraint (q 6 V  Þ, the optimal designs change more to closed wall cells as seen also in Table 1 [39]. The variety of designs presented in Table 1 covers a wide range of degrees of anisotropy as one can see through either the anisotropy plots or parameter A (defined in Section 2.3). From here one selects six designs to carry out the scale-size effects analysis: two with volume constraint only (shear-2 and shear-3) and

four with permeability constraint only (triaxial, biaxial, uniaxial and shear-3). Table 2 gives the corresponding homogenized stiffness tensors EH where coefficients are written in the order presented in Eq. (13) and Esolid =Ev oid ¼ 1012 is assumed. Tensors for triaxial, biaxial and uniaxial cases are clearly orthotropic while the other ones exhibit significant non-zero shear coefficients out of diagonal. Tensors for shear-3 cases are actually full tensors (21 non-zero coefficients) and include also negative coefficients (many  in the case of k ¼ 50%Þ. The numerical experiments, described in Sections 2.4, are carried out producing stiffness tensor estimates E based on global averages that are compared with the homogenized tensors EH . The bar charts in Figs. 6–11 summarize this comparative study. The ordinate axis is the percentage of deviation of an estimated Modulus, E , with respect to the homogenized value, EH , given by

Please cite this article in press as: Coelho PG et al. Scale-size effects analysis of optimal periodic material microstructures designed by the inverse homogenization method. Comput Struct (2015), http://dx.doi.org/10.1016/j.compstruc.2015.10.001

8

P.G. Coelho et al. / Computers and Structures xxx (2015) xxx–xxx 350

0 1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6

-35

E2323

E2233

Elastic coeficients

Elastic coeficients

(a)

(b)

50 25 0 -25

E2323

E3333

E2233

-50

E2323

75

E3333

100

1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6

E2233

125

600 550 500 450 400 350 300 250 200 150 100 50 0

E2222

150

Deviation [%]

1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6

175

E2222

E2323

0

E3333

-30

-75

1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6

-25

50

200

Deviation [%]

-20

E3333

100

-15

E2233

150

-10

E2222

200

-5

Deviation [%]

250

E2222

Deviation [%]

300

Elastic coeficients

Elastic coeficients

(c)

(d)

Fig. 7. Scale-size effects analysis for the biaxial case (k ¼ 50%). (a) Dirichlet-type BC’s with Esolid =Ev oid ¼ 1012 ; (b) Neumann-type BC’s with Esolid =Ev oid ¼ 10 ; (c) Constant stress on solid with Esolid =Ev oid ¼ 103 ; (d) Mixed BC’s with Esolid =Ev oid ¼ 103 . 

-20 -30 -40

Elastic coeficients

E3333

-50

E3333

0

Elastic coeficients

(a)

1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6

225 200 175 150 125 100 75

150 125 100 75 50

25

25

0

0

(c)

1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6

200 175

50

Elastic coeficients

(b)

250 225

E3333

5

-10

250

Deviation [%]

10

1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6

E3333

15

0

Deviation [%]

1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6

Deviation [%]

Deviation [%]

20

Elastic coeficients

(d)

Fig. 8. Scale-size effects analysis for the uniaxial case (k ¼ 50%). (a) Dirichlet-type BC’s with Esolid =Ev oid ¼ 1012 ; (b) Neumann-type BC’s with Esolid =Ev oid ¼ 10; (c) Constant stress on solid with Esolid =Ev oid ¼ 103 ; (d) Mixed BC’s with Esolid =Ev oid ¼ 103 . 

Eq. (18). The abscissa is the non-zero tensor component Eijkl and the different gray-colored bars refer to different scale factors used (n = 1 to 6).

Dev iation ½% ¼

  E  EH EH

 100

ð18Þ

The analysis of scale-size effects is focused on the porous materials thus the stiffness ratio Esolid =Ev oid ¼ 1012 is used and the homogenization values from Table 2 are used as targets. However, there are two exceptions where higher ratios are used in this study. Firstly, when alternative stress-based BC’s are applied (Fig. 4), one

uses Esolid =Ev oid ¼ 103 for designs with no connectivity in the solid phase along at least one direction (e.g. biaxial and uniaxial cases) in order to prevent rigid body motion of the solid phase. Secondly, one uses Esolid =Ev oid ¼ 10 for every experiment using Neumanntype BC’s in Eq. (7). Actually, this configures a composite mixing two solid materials rather than a porous one. In fact, Neumanntype BC’s are not effective for porous materials analyses due to the excessive compliance of the ‘‘void” phase near the boundary where the external pressure is applied. Analyses with Neumanntype BC’s are not skipped here, as seen in the charts below, but have to be carefully interpreted as referring to an ‘‘equivalent” bi-material composite solution. The comparative analysis with

Please cite this article in press as: Coelho PG et al. Scale-size effects analysis of optimal periodic material microstructures designed by the inverse homogenization method. Comput Struct (2015), http://dx.doi.org/10.1016/j.compstruc.2015.10.001

9

P.G. Coelho et al. / Computers and Structures xxx (2015) xxx–xxx 0

-20

Elastic coeficients

Elastic coeficients

(a)

(b)

E1313

E2323

E1212

E1223

E3313

E2213

E3333

E2222

E2233

E1113

E1313

E1223

E2323

E3313

E1212

E2213

E3333

E2222

E2233

E1113

E1122

E1133

-80

E1133

-60

1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6 E1111

-40

E1122

Deviation [%]

1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6

E1111

Deviation [%]

160 140 120 100 80 60 40 20 0 -20 -40 -60

Fig. 9. Scale-size effects analysis for the Shear-2 ðV  ¼ 50%Þ. (a) Dirichlet-type BC’s with Esolid =Ev oid ¼ 1012 . (b) Neumann-type BC’s with Esolid =Ev oid ¼ 10.

20

60

0

40

-20

Deviation [%]

0 -20 -40 -60 -80

1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6

-40 -60 -80 -100

1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6

-120 -140

-100 E1111 E1122 E1133 E1112 E1123 E1113 E2222 E2233 E2212 E2223 E2213 E3333 E3312 E3323 E3313 E1212 E1223 E1213 E2323 E2313 E1313

-160

E1111 E1122 E1133 E1112 E1123 E1113 E2222 E2233 E2212 E2223 E2213 E3333 E3312 E3323 E3313 E1212 E1223 E1213 E2323 E2313 E1313

Deviation [%]

20

Elastic coeficients

Elastic coeficients

(a)

(b)

125 1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6

1x1x1 2x2x2 3x3x3 4x4x4 5x5x5 6x6x6

100 75

Deviation [%]

50 25 0 -25 -50 -75 -100 -125 -150

E1111 E1122 E1133 E1112 E1123 E1113 E2222 E2233 E2212 E2223 E2213 E3333 E3312 E3323 E3313 E1212 E1223 E1213 E2323 E2313 E1313

400 350 300 250 200 150 100 50 0 -50 -100 -150 -200 -250

E1111 E1122 E1133 E1112 E1123 E1113 E2222 E2233 E2212 E2223 E2213 E3333 E3312 E3323 E3313 E1212 E1223 E1213 E2323 E2313 E1313

Deviation [%]

Fig. 10. Scale-size effects analysis for the Shear-3 case ðV  ¼ 50%Þ. (a) Dirichlet-type BC’s with Esolid =Ev oid ¼ 1012 . (b) Neumann-type BC’s with Esolid =Ev oid ¼ 10.

Elastic coeficients

Elastic coeficients

(a)

(b)

Fig. 11. Scale-size effects analysis for the Shear-3 case (k ¼ 50%). (a) Dirichlet-type BC’s with Esolid =Ev oid ¼ 1012 . (b) Neumann-type BC’s with Esolid =Ev oid ¼ 10. 

homogenization to be consistent requires updating of the tensors in Table 2 for different ratios, 101 or 103. The results plotted in Figs. 6–11 take already this update into account. Figs. 6–11 show the scale-size effects analysis based on the Dirichlet and Neumann BC’s applied to the designs selected from Table 2. Figs. 6–8 show also the results obtained with the alternative BC’s presented in Section 2.4.2, i.e. (1) Constant stress applied on solid phase only; (2) Uniform normal pressure upon rigid plate and shear based on displacement constraints (Mixed BC’s). Fig. 12 shows the convergence of the strain/stress energy densities evaluated from either Dirichlet or Neumann BC’s.

4. Discussion The results in Figs. 6–11 show that the diagonal tensor coefficients estimates applying Dirichlet-type BC’s converge always from above to the homogenization predictions as n increases. In contrast, the Neumann-type BC’s produce diagonal estimates that converge from below. This agrees with previous results reported in the literature [31–33]. However, this convergence trend may not be valid for the off-diagonal coefficients [31,32]. Although the orthotropic design cases treated here present the same trend for the diagonal and off-diagonal coefficients (see Figs. 6–8), the shear

Please cite this article in press as: Coelho PG et al. Scale-size effects analysis of optimal periodic material microstructures designed by the inverse homogenization method. Comput Struct (2015), http://dx.doi.org/10.1016/j.compstruc.2015.10.001

10

P.G. Coelho et al. / Computers and Structures xxx (2015) xxx–xxx

(a)

(b)

Fig. 12. Energy convergence for Dirichlet-type BC’s (dashed lines) and Neumann-type BC’s (solid lines): (a) strain energy density and (b) stress energy density.

design cases show a different trend (see Figs. 9–11). In fact, the inequalities in Eqs. (14) and (17) should be understood in the sense of tensors as defined in Eq. (15). The energy convergence charts in Fig. 12 confirm this. For the highest scale factor studied here (n = 6), all longitudinal coefficient estimates, i.e. EDiiii and ENiiii , are only 10% deviated from the homogenization results which is quite good for such a low factor. Also, for n = 6, the strain and stress energy densities evaluated from either Dirichlet or Neumann BC’s typically deviate less than  10% from the theoretical value (Shear-3 with k ¼ 50% and Dirichlet BC’s is here the exception). In [33] is shown poor convergence of Neumann estimates, ENiiii , because the application of the force to the compliant phase (Esolid =Ev oid ¼ 1012 was used) leads to high deformations in this phase near the boundary. This results in higher average strain and smaller elastic moduli predictions. In the present work, estimates of ENiiii target much better the homogenization values (deviation <10% for n = 6) because one has decided to replace here the porous composite by an ‘‘equivalent” bi-material composite with moderate contrast ratio, Esolid =Ev oid ¼ 10. The Neumann BC’s are then effective for non-porous composites. To treat porous composite materials, Neumann conditions can be slightly changed such that the uniform pressure is only applied on the solid phase. Having the analyses restricted to orthotropic cases (see Section 2.4.2), the resulting estimates deviate from homogenization less than 10%, as seen in Figs. 6c, 7c and 8c. The use of mixed BC’s improve these results even further, see Figs. 6d, 7d and 8d, most deviations are close to 0%. In general, the Dirichlet and Neumann BC’s seem to provide more conservative predictions while the alternative BC’s target much faster the theoretical values [33]. The present work addresses the rate of convergence of all elastic tensor coefficients as n increases. While the convergence of longitudinal coefficients is quite exceptional and consistent (regardless design and BC’s), that is not necessarily the pattern seen among the predictions for the other tensor coefficients. For example, in the case of orthotropic designs, non-longitudinal coefficients converge amazingly (deviations <7%) when Neumann conditions and the ratio Esolid =Ev oid ¼ 10 are used (see Figs. 6b, 7b and 8b). Furthermore, for this moderate contrast ratio one would get very pretty results as well using Dirichlet conditions, i.e. deviations <10% (also tested). However, when the porous case is treated (Esolid =Ev oid ¼ 1012 Þ, the rate of convergence is not so good for some

non-longitudinal coefficients. Even though, the alternative BC’s may provide better predictions comparing to Dirichlet’s (e.g. compare Fig. 6a with c and d, some deviations drop from 60% to 8%). Enlarging the number of targeted tensor coefficients up to 13 and 21 with the shear-2 and shear-3 cases, respectively, one sees actually an excellent convergence in Figs. 9a, b and 10a with most deviations less than 10% for n = 6. Larger deviations are obtained in Figs. 10b, 11a and b, especially in Fig. 11a. The Shear-3 case, with  k ¼ 50% and Dirichlet BC’s, gives the largest deviations in this study (can be up to 50% for n = 6) which is also evident in Fig. 12. From Table 2 one sees that the corresponding elastic tensor has comparatively greater complexity. Despite this, the estimate tensor ED agrees with EH in all coefficient signs and the highest deviations occur for coefficients that differ from the strongest moduli by one or two orders of magnitude. This difficult case treated here may show the importance of pursuing dedicated scale-size effects analysis for arbitrary anisotropic designs having n further increased. However, for most cases treated here, Fig. 12 shows a rapid rate of convergence of the apparent tensors to the corresponding homogenized tensors. Estimates of individual tensor coefficients (mainly non-longitudinal) may show greater variation in their convergence to the homogenized values. Finally, instead of averaging stresses and strains over all composite domain W one might averaging through the center unitcell domain Y only and then (as verified for all designs here) the apparent E and homogenization EH tensors would basically overlap (very small deviations, mostly up to 2% for n = 6). This means that the homogenization hypothesis come true in a local sense. Although one cannot exhaust all of the possible microstructure designs by studying only the six selected examples, one believes that this set can be quite representative as it is characterized by different anisotropy degrees: 0.38, 0.42, 0.46, 0.73, 1.08 and 2 (values of A in Table 1). In general, the results shown in Section 3 may suggest thus that for practical purposes it is enough to have a low scale factor (n = 6) to replace the non-homogeneous composite by the equivalent homogeneous material with the moduli predicted by homogenization theory.

5. Conclusions This work extends previous contributions on scale-size effects analysis of periodic material microstructures optimized by the

Please cite this article in press as: Coelho PG et al. Scale-size effects analysis of optimal periodic material microstructures designed by the inverse homogenization method. Comput Struct (2015), http://dx.doi.org/10.1016/j.compstruc.2015.10.001

P.G. Coelho et al. / Computers and Structures xxx (2015) xxx–xxx

inverse homogenization method (see e.g. [31–33]). The major contribution is the prediction of the full elastic tensor of threedimensional anisotropic designs confronting the homogenization result with the elastic properties assessed through standard testing procedures. The homogenization assumes periodic boundary conditions as well as infinite periodicity. In engineering practical materials and structures, this is not rigorously verified. Therefore, here we want to answer the question of how good homogenization predictions are compared to the mechanical responses of the actual composite obtained by numerical experiments. Firstly, one gets a set of optimal microstructures (unit-cells) solving a topology optimization problem that minimizes compliance subjected to permeability or volume constraints. Since one cannot exhaust all of the possible designs, one thinks that the few obtained designs here constitute a nice representative set for the scale-size effects analysis because it covers a large spectrum of anisotropy degrees. For the homogenization analysis, only one unit-cell is analyzed since the periodicity assumption guarantees the same result for any number of cells. The actual composite is generated by an array of n  n  n unit-cells, here n = 1 to 6. The resulting composite for each scale factor ‘‘n” is subjected to six basic numerical tests (three normal and three shear tests) applying stress or displacementbased BC’s. The mechanical responses are measured in terms of stresses and strains. Ratios between averages provide elastic coefficient estimates which are n-dependent, and the objective is to investigate their convergence to the homogenization values. In general, porosity brings convergence challenges especially when applying Neumann-type conditions due to the excessive compliance of the void phase where a constant pressure is applied. When one replaces the void phase by a solid one (one gets a bimaterial composite), the convergence improves significantly. One also observes that faster convergence is obtained for longitudinal coefficients and when alternative BC’s are used instead of standard Neumann and Dirichlet BC’s. However, these alternative BC’s are only applied for the orthotropic design cases where static equilibrium can be easier achieved. Among all the designs here, the offdiagonal and non-longitudinal coefficients show a slower convergence and this is evident in one out of six designs studied. This may justify careful and auxiliary scale-size effects analysis for other arbitrary anisotropic designs. In general, strain and stress energy measures show faster convergence than individual tensor coefficients. In most cases, resulting errors drop to 5 to 10% for a small number of unit-cell repetitions. The present work indicates that for practical purposes, even for micro-structures with low scale factor (n = 5 or 6), it is mechanically admissible to use the equivalent moduli given by the periodic homogenization theory. These observations are also consistent with previous works focusing on bi-dimensional microstructures with material symmetry. Acknowledgments This work was partially supported by FCT, through IDMEC, under LAETA, project UID/EMS/50022/2013FCT-Portugal and through UNIDEMI, project UID/EMS/00667/2013FCT-Portugal. References [1] Cadman JE, Zhou S, Chen Y, Li Q. On design of multi-functional microstructural materials. J Mater Sci 2013;48:51–66. [2] Anfeassen E, Lazarov BS, Sigmund O. Design of manufacturable 3D extremal elastic microstructure. Mech Mater 2014;69:1–10. [3] Sigmund O. Systematic design of metamaterials by topology optimization. In: IUTAM symposium on modelling nanomaterials and nanosystems. vol. 13; 2009. p. 151–9. [4] Sigmund O. Materials with prescribed constitutive parameters: an inverse homogenization problem. Int J Solids Struct 1994;31(17):2313–29.

11

[5] Bendsøe M, Sigmund O. Topology optimizarion. theory, methods and applications. Berlin: Springer; 2003. [6] Andreasen CS, Sigmund O. Multiscale modeling and topology optimization of poroelastic actuators. Smart Mater Struct 2012;21:065005. [7] Huang X, Zhou SW, Xie YM, Li Q. Topology optimization of microstructures of cellular materials and composites for macrostructures. Comput Mater Sci 2013;67:397–407. [8] Yang XY, Huang X, Rong JH, Xie YM. Design of 3D orthotropic materials with prescribed ratios for effective Young’s moduli. Comput Mater Sci 2013;67:229–37. [9] Radman A, Huang X, Xie YM. Topological optimization for the design of microstructures of isotropic cellular materials. Eng Optim 2013;45 (11):1331–48. [10] Radman A, Huang X, Xie YM. Topological design of microstructures of multiphase materials for maximum stiffness or thermal conductivity. Comput Mater Sci 2014;91:266–73. [11] Challis VJ, Guest JK, Grotowski JF, Roberts AP. Computationally generated cross-property bounds for stiffness and fluid permeability using topology optimization. Int J Solids Struct 2012;49:3397–408. [12] Zuo ZH, Huang X, Rong JH, Xie YM. Multiscale design of composite materials and structures for maximum natural frequencies. Mater Des 2013;51:1023–34. [13] Xie YM, Yang X, Shen J, Yan X, Ghaedizadeh A, Rong J, et al. Designing orthotropic materials for negative or zero compressibility. Int J Solids Struct 2014;51:4038–51. [14] Andreasen CS, Andreassen E, Jensen JS, Sigmund O. On the realization of the bulk modulus bounds for two-phase viscoelastic composites. J Mech Phys Solids 2014;63:228–41. [15] Chen W, Liu S. Topology optimization of microstructures of viscoelastic damping materials for a prescribed shear modulus. Struct Multidisc Optim 2014;50:287–96. [16] Andreassen E, Jensen JS. Topology optimization of periodic microstructures for enhanced dynamic properties of viscoelastic composite materials. Struct Multidisc Opim 2014;49:695–705. [17] Dong-ming X, Yong-qiang Y, Xu-bin S, Di W, Zi-yi L. Topology optimization of microstructure and selective laser melting fabrication for metallic biomaterial scaffolds. Trans Nonferrous Met Soc China 2012;22:2554–61. [18] Vatanabe SL, Paulino GH, Silva ECN. Design of functionally graded piezocomposites using topology optimization and homogenization – toward affective energy harvesting materials. Comput Meth Appl Mech Eng 2013;266:205–18. [19] Jayachandran KP, Guedes JM, Rodrigues HC. Piezoelectricity enhancement in ferroelectric ceramics due to orientation. Appl Phys Lett 2008;92(23):232901. [20] Bourgat JF. Numerical experiments of the homogenization method for operators with periodic coefficients. Lect Notes Math 1977;704:330–56. [21] Bensoussan A, Lions JL, Papanicolaou G. Asymptotic analysis for periodic structures. North-Holland; 1978. [22] Sanchez-Palencia E. Lect Notes Phys, vol. 127. Springer-Verlag; 1980. [23] Guedes JM, Kikuchi N. Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element method. Comput Meth Appl Mech Eng 1990;83:143–98. [24] Zohdi TI, Wriggers P. Introduction to computational micromechanics. Springer; 2004. [25] Zuo ZH, Huang X, Yang X, Rong JH, Xie YM. Comparing optimal material microstructures with optimal periodic structures. Comput Mater Sci 2013;69:137–47. [26] Xie YM, Zuo ZH, Huang X, Rong JH. Convergence of topological patterns of optimal periodic structures under multiple scales. Struct Multidisc Optim 2012;46:41–50. [27] GaoMing D, WeiHong Z. Size effects of effective Young’s modulus for periodic cellular materials. Sci China Ser G-Phys Mech Astron 2009;52(8):1262–70. [28] Dai G, Zhang W. Cell size effects for vibration analysis and design of sandwich beams. Acta Mech Sin 2009;25:353–65. [29] Huang X, Xie YM. Optimal design of periodic structures using evolutionary topology optimization. Struct Multidisc Optim 2008;36:597–606. [30] Zhang W, Sun S. Scale-related topology optimization of cellular materials and structures. Int J Numer Meth Eng 2006;68:993–1011. [31] Pecullan S, Gibiansky LV, Torquato S. Scale effects on the elastic behavior of periodic and hierarchical two-dimensional composites. J Mech Phys Solids 1999;47:1509–42. [32] Hollister SJ, Kikuchi N. A comparison of homogenization and standard mechanics analyses for periodic porous composites. Comput Mech 1992;10:73–95. [33] Coelho PG, Guedes JM, Rodrigues HC. Scale-size effects analysis on predicting the mechanical properties of periodic composite materials. In: Topping BHV, Iványi P, editors. Proceedings of the twelfth international conference on computational structures technology. Stirlingshire, UK: Civil-Comp Press; 2014. http://dx.doi.org/10.4203/ccp.106.247. Paper 247. [34] Dias MR, Guedes JM, Flanagan CL, Hollister SJ, Fernandes PR. Optimization of scaffold design for bone tissue engineering: a computational and experimental study. Med Eng Phys 2014;36(4):448–57. [35] Coelho PG, Fernandes PR, Rodrigues HC. Multiscale modeling of bone tissue with surface and permeability control. J Biomech 2011;44(2):321–9. [36] Kang H, Lin C-Y, Hollister SJ. Topology optimization of three dimensional tissue engineering scaffold architectures for prescribed bulk modulus and diffusivity. Struct Multidisc Optim 2010;42:633–44.

Please cite this article in press as: Coelho PG et al. Scale-size effects analysis of optimal periodic material microstructures designed by the inverse homogenization method. Comput Struct (2015), http://dx.doi.org/10.1016/j.compstruc.2015.10.001

12

P.G. Coelho et al. / Computers and Structures xxx (2015) xxx–xxx

[37] Xu S, Cheng G. Optimum material design of minimum structural compliance under seepage constraint. Struct Multidisc Optim 2010;41:575–87. [38] Guest J, Prévost JH. Optimizing multifunctional materials: design of microstructures for maximized stiffness and fluid permeability. Int J Solids Struct 2006;43:7028–47. [39] Sigmund O. On the optimality of bone microstructure. In: Pedersen P, Bendsøe MP, editors. IUTAM symposium on synthesis in bio solid mechanics. Kluwer Academic Publishers; 1999. p. 221–34. [40] Guedes JM, Rodrigues HC, Bendsøe MP. A material optimization model to approximate energy bounds for cellular materials under multiload conditions. Struct Multidisc Optim 2003;25:446–52. [41] Coelho PG, Rodrigues HC. Hierarchical topology optimization addressing material design constraints and application to sandwich-type structures. Struct Multidisc Optim 2015. http://dx.doi.org/10.1007/s00158-014-1220-x.

[42] Bendsøe MP. Optimal shape design as a material distribution problem. Struct Optim 1989;1:193–202. [43] Svanberg K. The method of moving asymptotes – a new method for structural optimization. Int J Numer Meth Eng 1987;24:359–73. [44] Hashin Z, Shtrikman S. A variational approach to the theory of effective magnetic permeability of multiphase materials. J Appl Phys 1962;33:3125–31. [45] Milton GW, Kohn RV. Variational bounds on the effective moduli of anisotropic composites. J Mech Phys Solids 1988;36:597–629. [46] Böhlke T, Brüggemann C. Graphical representation of the generalized Hooke’s Law. Tech Mech 2001;21:145–58. [47] Challis VJ, Roberts AP, Wilkins AH. Design of three dimensional isotropic microstructures for maximezed stiffness and conductivity. Int J Solids Struct 2008;45:4130–46.

Please cite this article in press as: Coelho PG et al. Scale-size effects analysis of optimal periodic material microstructures designed by the inverse homogenization method. Comput Struct (2015), http://dx.doi.org/10.1016/j.compstruc.2015.10.001