Stress localization and strength optimization of frame material with periodic microstructure

Comput. Methods Appl. Mech. Engrg. 197 (2008) 4016–4026

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

Stress localization and strength optimization of frame material with periodic microstructure Fabian Lipperman *, Moshe B. Fuchs, Michael Ryvkin School of Mechanical Engineering, Tel Aviv University, Tel Aviv 69978, Israel

a r t i c l e

i n f o

Article history: Received 3 December 2007 Received in revised form 25 February 2008 Accepted 26 March 2008 Available online 3 April 2008 Keywords: Maximum strength Cellular materials Localization Homogenization Hashin–Shtrikman bounds

a b s t r a c t Stress localization for an elastic material with periodic microstructure is applied to the optimal design of rigid-jointed cellular material. The discretized form of localization was published some 20 years ago but was primarily used for strain localization (and stiffness homogenization). The stress localization aspect did not receive much attention in the past and even more so for materials modeled as rigid- or pinjointed lattices. With the present minimization of the maximum local stress under a macroscopic stress state efficient stress localization was a major concern and it has therefore received due consideration in this paper. A ground-structure approach was used for the optimal design where in addition to topological variables for the unit cell geometric variables determine the coordinates of two internal nodes and the overall shape of the representative volume element. In a first instance the material was subjected to a uniaxial and a pure shear macroscopic state of stress of arbitrary orientation. Three regular shapes emerged from the optimization: a square crossed by diagonals, an equilateral triangle and the Kagomé. The three figures, which exhibit predominant axial modes with no significant bending, have very similar strengths. In a second example the ground-structure was subjected to a transversely isotropic macroscopic tensile stress and to a transversely isotropic macroscopic compressive stress. The design space had a plethora of optimal configurations of same maximum strength. It is shown that these optimal designs, although anisotropic, expand isotropically under isotropic stresses and reach the Hashin–Shtrikman upper bounds for the bulk modulus of isotropic materials in the low density range. In the compressive case, slender elements were eliminated by the Euler buckling constraint and bulkier optimal patterns were obtained. Ó 2008 Elsevier B.V. All rights reserved.

1. Introduction The subject matter of this paper is the design for maximum strength of a material with periodic microstructure that is modeled as an infinite repetitive structure composed of rigid-jointed uniform elastic elements. Owing to the periodicity of the model the redesign can be confined to a single repetitive cell and one looks for a distribution of material between the uniform elements which will minimize the maximum tensile stress anywhere in the cell when subjected to one or more macroscopic stress fields. In view of the multiple structural reanalysis implied in any numerical optimization algorithm, an efficient method for analyzing the cellular material when in a macroscopic stress field is a major concern. We recall that when an infinite repetitive structure is in a uniform macroscopic state of stress it can be assumed that the local (or detailed) stresses are periodic. Obtaining the local stresses in a unit cell on the basis of the macroscopic stresses is called stress localization. (Here and in the sequel we follow the nomenclature * Corresponding author. Tel.: +972 36405579; fax: +972 36407617. E-mail address: [email protected] (F. Lipperman). 0045-7825/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2008.03.019

in Suquet [1].) Similarly, strain localization is the procedure for determining the local strains in the unit cell when the material is subjected to a uniform macroscopic strain field. Arguably the most important implementation of localization consists in obtaining the effective or macroscopic properties of the material, that is, the relationship between the macroscopic stresses and macroscopic strains. This is termed homogenization and indeed, for a twodimensional (2D) material, localizations for three unit macroscopic strain fields lead to the effective or macroscopic elastic stiffnesses of the material. Similarly, three stress localizations will produce the macroscopic elastic compliances. Numerical examples of strain localizations are legion, so much so, that the term ‘homogenization’ became akin to ‘stiffness homogenization’. Interestingly, stress localization appeared as the alter-ego of strain localization in conference proceedings as early as 1985 [2] and in later publications by the same research group [3] however, applications, although extant [4], are relatively scarce. This state of affairs prevails for the linearly elastic materials with framed microstructure, considered herein. These are 2D cellular materials modeled either with bars in pin-jointed (truss) or with beams in rigid-jointed (frame) models (see [5]). Strain localization

4017

F. Lipperman et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4016–4026

with periodic boundary conditions was used by Sigmund [6,7] for designing the topology of frame materials with optimal macroscopic elastic stiffnesses. For other applications based on strain localization see [8–13] for recent publications. All followed a similar path for computing the macro-stiffnesses of frame or truss material. We should also mention [14] where macroscopic stiffnesses were obtained as a by-product of the study of macroscopic mechanisms of periodic pin-jointed trusses. Unlike strain localization for which standard techniques exist, stress localization seems, in practice, to be an ad hoc affair. In the particular cases of cellular materials with relatively simple geometries, for instance lattices with isotropic or square symmetry, one can infer the exact stress conditions at the cell boundaries (see, for instance, [15–18]). However, in the absence of at least two axes of symmetry, the exact boundary conditions cannot always be defined and stress localization becomes problematic. As noted (see [9,13,19]), in the case of grid structures with high structural complexity, for exact results, [stiffness] homogenization with periodic boundary conditions is the way to go. If strain localization is, more often than not, used as a tool for computing the macroscopic elastic stiffnesses, stress localization can be a goal in itself. Indeed, for determining the strength (for instance, the maximum local stress, as done herein) of a cellular material in a macroscopic stress state, there is no need for the effective properties, which require three (usually strain) localizations, that is, three analyses. A single stress localization (one analysis) will provide the answer. For the analysis we will therefore use the discretized form of localization given in [3] which is very general and applies equally to stress and strain localization and for unit cells with or without internal symmetry. With regard to the optimal design of cellular materials we note that previous publications deal primarily with optimizing macroscopic (effective) quantities [6,7,13,20,21]. As indicated, in this work we design cellular materials for maximum strength although the maximum effective bulk modulus will appear, as an afterthought, in one of the numerical examples. The material is modeled as an infinite repetitive structure composed of rigid-jointed uniform elastic members and the strength of the material, in the fashion of the weakest link, is given by the maximum stress anywhere in the unit cell, either rupture or buckling. A modified ground-structure approach [7,22–24] is used for the optimization where from an initial grid of redundant structures the optimal topology of the repetitive cell emerges by elimination of ineffective members and reshaping of the cell boundaries. The design variables are the cross-sectional areas of the elements of the unit cell augmented by design parameters controlling the shape of the cell and the coordinates of two internal nodes. Two cases were investigated. In a first instance, the material is in a state of macroscopic uniaxial stress and in a state of macroscopic pure shear. The direction of the loading is arbitrary and this multiple load conditions problem was solved by considering in sequence a dense set of orientations covering half a circle. In the second example the material was subjected to a macroscopic transversely isotropic tensile stress state and to a macroscopic transversely isotropic compressive state of stress. The outline of the paper is as follows. Stress and strain localization is presented in Section 2. For want of repeating a very clear discourse published in [3] we will describe an apparently shorter, albeit less rigorous, route for reaching the discretized form of localization. In Section 3 the technique is specialized for cellular materials in the form of infinite grids and as an example localization is applied to a rigid-jointed model of the triangle–triangle grid (T–T), a structure analyzed by Hutchinson and Fleck [14] for a pin-jointed model. In the subsequent section two optimal design problems are introduced and the results are described and discussed in Section 5. The paper ends with a summary in the last section.

2. Localization with periodic boundary conditions We present in this section a general displacement-based technique of localization for linear-elastic materials with arbitrary repetitive micro-structures, subjected to either a macroscopic stress field (stress localization) or a macroscopic strain field (strain localization). As noted, the theory has been available for many years [2,3] but we will confine ourselves to presenting a slightly different route to the localization equations with the necessary particularization for rigid-jointed grids. The representative volume element (RVE), a single unit cell with periodic boundary conditions, is composed of beam elements. Without loss of generality, this section will first consider rectangular unit cells. Subsequently, the method will be generalized to parallelogram unit cells. Consider a 2D material, formed by the repetition of a rectangular unit cell (Fig. 1a), in a macroscopic stress/strain field r ¼ f rx

ry

sxy gT ;

e ¼ f ex

ey

cxy gT ;

ð1Þ

where the components of the macroscopic stresses and strains are given with respect to orthogonal axes x and y. The generic domain is that part of the RVE which by translations can generate the infinite structure. It will be noted that the entire RVE except the right (V+) and upper (H+) borders, including their end-nodes, is the generic domain. The nodes of the RVE, schematically depicted in Fig. 1a, are correspondingly divided into generic nodes (black dots) and non-generic ones (white dots). The left (V) and bottom (H) borders are therefore generic (in continuous lines). We start by assembling the stiffness matrix K0 of the RVE (line elements along the non-generic borders are not assembled). This leads to equilibrium equations of the form K0 u0 ¼ p0 ;

ð2Þ

where u0 is the displacements of all the nodes of the RVE and p0 contains the corresponding nodal loads. For convenience we assume that the nodal sequence in u0 and correspondingly in p0 is u0 ¼ f u0I

u0H

u0V

u0H V

u0Hþ

u0Vþ

u0H Vþ

u0Hþ V

u0Hþ Vþ gT ;

f p0I

p0H

p0V

p0H V

p0Hþ

p0Vþ

p0H Vþ

p0Hþ V

p0Hþ Vþ

ð3Þ 0

p ¼

T

g ; ð4Þ

u0I u0Hþ

where is the displacements of the internal nodes, u0H u0V u0Vþ are the displacements of the nodes on the borders of the RVE, excluding the corner nodes, and u0H V u0H Vþ u0Hþ V u0Hþ Vþ are the four corner nodes displacements. The corner nodes vectors have each three components whereas the lengths of the other displacement vectors along the boundaries are multiples of three depending on the number of nodes on the borders. One will note that the first four entries of u0 relate to the generic nodes while the non-generic ones are to be found in the last five components of u0 . The force vector p0 is partitioned accordingly. To the equilibrium equations (2) one must add the periodicity conditions on the boundaries of the RVE [1]. For the displacement periodicity conditions we note that since the solution is periodic the deformations of all unit cells are identical. In 2D, that deformation can be represented by (see Fig. 2) 2 3 LH 0 0 6 7 d ¼ Le with L ¼ 4 0 LV 0 5 and d ¼ f dx dy dxy gT ; ð5Þ 0

0

LV

where LH and LV are the base and height of the unit cell, respectively, dx and dy are the change in length of the horizontal and vertical boundaries, respectively, and dxy is the shift of the top horizontal boundary with respect to the bottom one (Fig. 2b) Note, cxy is the engineering strain (cxy = 2exy). It is easy to see that with the

4018

F. Lipperman et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4016–4026

H+V-

H+ V-

H+ V+

H+

H+V+

H+

n θn

y V-

V+

-

HV

H-

-

-

V-

x

V+

HV

H- V+

H-

H-V-

+

a

b

Fig. 1. (a) A rectangular RVE and (b) a parallelogram-shaped RVE.

δxy LH

LH + δx c'

d'

d

d

c

c LV + δy

y

LV

x a

b

a, a'

b b'

a

b

Fig. 2. The modes of deformation of the rectangular unit cell: (a) an original unit cell before; (b) after deformation.

H+ dyj+δy dxj+δxy

H+V-

the node. With reference to (2), we find that the nine components of u0 are related to the five components of a generalized displacements vector (Fig. 3)

H+V+

θj

u ¼ f uI +

V

Vdyi dxi

θi -

HV

dyi θi

dyj θj

-

dxj

dxi+δx

-

HV

+

-

H

Fig. 3. The degrees of freedom of generic and non-generic nodes in an RVE.

displacements periodicity conditions the displacements of repeating nodes in adjacent unit cells are augmented by the components of d [25], as shown in Fig. 3. For instance, the degrees of freedom of the nodes along V+ and H+, u0iVþ and u0jHþ , will be related to the corresponding nodes along V and H by u0iVþ ¼ u0iV þ f dx u0jHþ

¼

u0jH

þ f dxy

0 dy

0 gT ; 0 gT

ð6Þ

assuming that the displacements of a typical node i are in the sequence u0i ¼ f dix

i

dy

T hi g ;

ð7Þ

where dx and dy are translational displacements in the x and ydirections, respectively, and h is the counterclockwise rotation of

uH

uV

by 9 2 8 u0I > I > > > > 6 > > > > u0H > 0 6 > > > > > 6 > > 60 > > u0V > > > 6 > > > > > > 6 0 > > 0 > = 6 < uH V > 6 0 uHþ ¼6 0 6 > > > 6 > > u0Vþ > > > 60 > > > > > 6 > 0 > 60 uH Vþ > > > > > 6 > > > > 6 0 > > > uHþ V > 40 > > > > ; : 0 uHþ Vþ 0

0 I

uH V

d gT

0 0

0 0

0 0

0

I

0

0

0

0

I

0

I

0

0

BH

0

I

0

BV

0

0

I

BV

0

0

I

BH

0

0

I

BH þ BV

ð8Þ 3 7 7 9 78 7> uI > 7> > > > > 7> > > uH > 7> = 7< 7 7> uV > > 7> > > uH V > 7> > > > 7> 7: d ; 7 7 5

ð9Þ

or symbolically u0 ¼ Gu:

ð10Þ

Eq. (10) is the displacement periodicity conditions for the discretized model. The dimensions of the zero (0) and unit matrices (I) in G are commensurate with their position in the matrix. Likewise, matrices BH and BV are formed of submatrices 2 3 2 3 0 0 1 1 0 0 6 7 6 7 BH ð0Þ ¼ 4 0 1 0 5 and BV ð0Þ ¼ 4 0 0 0 5; ð11Þ 0

0

0

0

0

0

respectively, stacked one on top of the other, the number of submatrices in any BH and BV depending on their location in matrix G. Using the change of variables in (10) produces the final form of the equilibrium equations

4019

F. Lipperman et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4016–4026

Ku ¼ p;

"

ð12Þ

where the reduced stiffness matrix is given by the congruent transformation K ¼ GT K0 G and ð14Þ

The congruent load vector p can be partitioned in sub-vectors which correspond to the sub-vectors in u (8) p ¼ f pI

pH

pV

pH V

u gT :

Kgd

KTgd

Kdd

ug d



 ¼

0



ðaÞ; ðbÞ;

u

ð19Þ

where ug is the displacements of the generic nodes. Eq. (19) governs the analysis (localization) of the RVE. (It is equivalent to Eq. (45) in [3].) When employed for strain localization, that is, when d is given, Eq. (19a) is solved for ug. The generic displacements in conjunction with d give the displacements of all the other nodes of the RVE (10). When stress localization is the purpose (u is given) system (19) is solved and with ug and d at hand one proceeds as for the strain localization. Before concluding this section it is perhaps instructive to remind the reader that the above procedure rests on the assembly of two standard matrices K0 and G. The former is the classical stiffness matrix of the generic elements of the RVE and the latter is a simple matrix expressing the displacements of the nodes of the RVE in terms of the displacements of the generic nodes. Noteworthy, the present method does not explicitly divide the strains in average and fluctuating terms as is customary in homogenization theory (see [3]). It, arguably, simplifies the numerical implementation of localization, especially when assembling G (10).

ð13Þ

p ¼ GT p0 :

Kgg

#

ð15Þ

It is instructive to expand (14) 9 9 8 8 p0I pI > > > > > > > > > > > > > > > > 0 0 > > > p  þp þ  > p > > > > H H H = = < < 0 0 pV þ pVþ : ¼ pV > > > > > > > > > > p0H V þ p0H Vþ þ p0Hþ V þ p0Hþ Vþ > > > > > pH V > > > > > > > > ; ; > : : BV ðp0Vþ þ p0H Vþ þ p0Hþ Vþ Þ þ BH ðp0Hþ þ p0Hþ V þ p0Hþ Vþ Þ u ð16Þ We now introduce the stress periodicity conditions which stipulate that the tractions on opposite borders of the RVE are equal and opposite. For the discretized model we get

3. Particularization to reticulated structures 3.1. Parallelogram-shaped RVE

p0Hþ ¼ p0H ; p0Vþ ¼ p0V ; p0Hþ Vþ ¼ p0H V ;

Rectangular RVEs with their orthogonal axis seem an easy way to convey the basic ideas underlying localization for a discretized model. However, unit cells in the form of parallelograms are often used [26] and are a natural choice in many cases (see examples in [14,25]). In Fig. 4 we show parallelogram-shaped unit cells and their generic and non-generic nodes for three common cellular materials with hexagonal, triangular and Kagomé layouts. An RVE in the form of a parallelogram was also shown in Fig. 1b. As in the rectangular cases, H and H+ are the lower and upper boundaries of the RVE and V and V+ its left and right lateral boundaries, respectively. We denote by n the direction of the outer normal to the inclined border V+. The slope is defined by hn which is the counterclockwise rotation of axis n relative to the x-axis. The deformation of the unit cell is here conveniently given by

ð17Þ

p0Hþ V ¼ p0H Vþ ; which when introduced in (16) leads to pH ¼ pV ¼ pH V ¼ 0. Now, pI = 0 because we assume that there are no body forces and, consequently, the load vector is zero but for u. The latter is forces congruent with the deformations d of the RVE. A perusal of the last line in (16) will evidence that the components of u (ux, uy, uxy) are simply the sum of all the forces applied to border V+ of the RVE in the x-direction, the sum of all the forces applied to border H+ in the y-direction and the sum of all the forces applied to border H+ in the x-direction, respectively. These can be given in terms of the macroscopic stress field: 2 3 LV 0 0 6 7   u ¼ L r with L ¼ 4 0 LH 0 5: ð18Þ 0 0 LH

d ¼ f dn

dxy gT ;

ð20Þ

where dn, the deformation of the unit cell along the n direction, replaces dx of the rectangular case. One assembles the equilibrium equations (2) in a standard manner, however, for the reduction of the degrees of freedom via (10) submatrix BV(0) is replaced by

The reduced equilibrium equations can conveniently be partitioned as

a

dy

b

c

Fig. 4. Examples of parallelogram-shaped unit cells composing the (a) hexagonal, (b) triangular, and (c) Kagomé lattices. Black nodes are generic, whites are non-generic.

4020

F. Lipperman et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4016–4026

2

cos hn 6 BV ðhn Þ ¼ 4 sin hn 0

0

0

3

0

7 0 5:

0

0

yielding equilibrium equations of the form 2 38 9 8 9 > > >0> > > > u2 > > > > > > > > 6 7> > 0> u3 > > > > > 6 7> = < = < 6 7 6 7 ð25Þ K 6 7> u4 > ¼ > 0 > > > > 6 7> > 0> u1 > > > > > 4 5> > > > > > > ; > : > ; : u d pffiffiffi with u ¼ ry 3Lf1=4; 1; 0gT . In the above equations the zero and unit matrices are 3  3 and the zero vectors are 3  1. To prevent rigid body motion we set the two translational degrees of a generic node to zero. In view of the periodicity the rigid-body rotation is thus also zero. This leads to a linear system of 15  2 = 13 unknown degrees of freedom from which the stress resultants (bending moments, shear and axial forces) of the elements of the RVE are readily computed. The results of the stress localization using Bernoulli–Euler beams are given in Table 1, with g = h/L, where node a is the lower node of each element (see Fig. 5b). To conclude this example we have also computed the macroscopic elastic compliances D by performing similar localizations for rx = 1 and sxy = 1 (in addition to ry = 1). We solve the reduced Eq. (25) with the corresponding components of u : u1 ¼ oT pffiffiffi pffiffiffi npffiffiffi 3Lf3=4; 0; 0gT and u3 ¼ 3L 3=2; 0; 1 . This leads to 9 9 8 0 18 pffiffiffi rx > D1 D2 0 > > = = < ex > < 3 B C ey ¼ @ D2 D1 0 A ry > > > > 12ED 0 ; ; : : cxy sxy 0 0 D3 9 9 8 8 0 1 C1 C2 0 > > = pffiffiffi = < rx > < ex > 3Eg B C ry ¼ ð26Þ @ C 2 C 1 0 A ey > > > > 2C 0 ; ; : : cxy sxy 0 0 C3

ð21Þ

Going through the algebra we obtain the reduced equilibrium equations (19) where the components of u are now 8 9 > < rn LV > = u ¼ ry LH ; > > : ; sxy LH

ð22Þ

with rn ¼

rx þ ry rx  ry þ cos 2hn þ sxy sin 2hn : 2 2

ð23Þ

Note, a direct assembly of the reduced stiffness matrix K, without passing through K0 and G, is also possible. 3.2. Example of stress localization As an example of localization we will analyze the triangular– triangular (T–T) infinite structure subjected to a macroscopic uniaxial tensional stress in the y-direction, ry(rx = sxy = 0), shown in Fig. 5a. This cellular material was analyzed for a pin-jointed model in [14] and we will give the rigid-jointed solution using stress localization. The selected RVE (Fig. 5b) encloses four generic nodes (black dots), three non-generic nodes (white dots), and seven uniform elements with widths h and unit inward depth. The elements are of length L but for the last two which have length L/2. (Elements 6 and 7 are in fact halves of a same element which crosses the boundaries). The slopes of the elements with respect to the x-axis are multiples of p/6 and Young’s modulus of the solid material is E. Noting that for the T–T pattern u0I ¼ fu2 u3 gT ; u0H ¼ u4 ; u0V ¼ hvoidi; u0H V ¼ u1 ; u0Hþ ¼ u6 ; u0Vþ ¼ hvoidi; u0H Vþ ¼ hvoidi; u0Hþ V ¼ u5 and u0Hþ Vþ ¼ u7 the transformation (10) becomes 8 9 2 I u2 > > > > > > 60 > > > > u 6 > > 3 > 6 > > > 6 > > > > > = 60 < u4 > 6 u1 ¼ 6 0 > > 6 > > > > 6 > > 60 > u6 > > > > > 6 > > > 40 > u5 > > > ; : > u7 0

with Table 1 The internal axial forces and end-moments of the elements of the T–T rigid-jointed pattern under a macroscopic uniaxial tensional stress ry (units are in parenthesis)  pffiffi    3r L ry L2 Element Axial force 4ðg2 þ2Þðgy 2 þ1Þ End moments 8ðg2 þ2Þðg 2 þ1Þ

3

0

0

0

0

I

0

0

0

0 0

I 0

0 I

0 0

0

I

0

BH ð0Þ

0

0

I

BH ð0Þ

0

0

I

BH ð0Þ þ BV ðp=6Þ

78 u 9 7> 2> > > 7> > > > > 7> u > = 7< 3 > 7 7 u4 > 7> > > u1 > 7> > > > 7> > ; 7: > 5 d

g4 + 4g2 + 2 g4 + 4g2 + 2 g4 + 6g2 + 6 g4  g2  2 g4  g2  2 g4 + 6g2 + 6 g4 + 6g2 + 6

1 2 3 4 5 6 7

ð24Þ

5

δxy 6 6

(6g4 + 9g2 + 2) (4g4 + 7g2 + 2) (3g2 + 2) (g4 + g2 + 2) (3g4 + 3g2  2) (g4 + g2) (2g4  g2  2)

7 4

5

L V = 3L

2 2

y

7

n θn = π / 6

4

1

x

δn

3

3 1

a

Node b

δy

σy

σy

Node a (4g4 + 7g2 + 2) (6g4 + 9g2 + 2) (2g4  g2  2) (3g4 + 3g2  2) (g4 + g2 + 2) (3g2 + 2) (g4 + g2)

L H = 3L

b

Fig. 5. (a) The triangular–triangular (T–T) infinite structure subjected to a macroscopic uniaxial tensional stress and (b) its repetitive parallelogram.

4021

F. Lipperman et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4016–4026

D0 ¼ g7 þ 3g5 þ 2g3 ; D1 ¼ 3g6 þ 24g4 þ 25g2 þ 2;

C 0 ¼ 21g4 þ 31g2 þ 8; C 1 ¼ 3g6 þ 24g4 þ 25g2 þ 2;

D2 ¼ 3g6  4g4  7g2 þ 2;

C 2 ¼ ð3g6  4g4  7g2 þ 2Þ;

D3 ¼ 56g4 þ 64g2 ;

C 3 ¼ 3g6 þ 10g4 þ 9g2 þ 2:

ð27Þ

We have added the macroscopic stiffnesses C to show that when g approaches zero the effective stiffnesses of (26) coincide with a relation obtained for trusses [14]: 9 8 0 1 > = pffiffiffi < rx > 3Eg B ry ¼ @ 1 > > 8 ; : sxy 0

1 1 0

nates of the internal nodes with respect to the center of the parallelogram. The area A of the parallelogram A ¼ LH LY ¼ LH LV cos hn

ð30Þ

is fixed and consequently the cell aspect ratio, LV/LH, and the angle of the lateral boundaries, hn, depend on the design variables LY and LV by LV =LH ¼ LV LY =A;

ð31Þ

hn ¼ cos1 ðLY =LV Þ:

9 18 0 > = < ex > C 0 A ey : > > ; : cxy 1

ð28Þ

4. Cellular materials of maximum strength

Recalling that the objective of the design is to reduce the highest axial stress anywhere in the RVE when the material is in one or more macroscopic stress states, we obtain the following min–max optimization problem min max rle

ðaÞ

v

With stress localization at hand we have an efficient method (one analysis) for determining the strength of cellular material for given states of macroscopic stress. Turning our attention to the design we will minimize the maximum axial stress anywhere in the unit cell for all considered states of macroscopic stress while preventing the occurrence of local buckling in compressive elements. We are using a ground-structure approach similar to what was done in [7] but instead of optimizing macroscopic effective properties we are minimizing local stresses. In other words, the analysis module is concerned with stress localization rather then homogenization. The design variables are the thicknesses of the uniform beam elements connecting the nodes of the RVE augmented by parameters which control the geometry of the RVE parallelogram (see [21]) and the coordinates of two internal nodes. We have thus a combination of topology and geometry variables to enhance the solution. Also, no symmetry of the microstructure was assumed. The RVE, shown in Fig. 6, has five generic nodes (1–5) and 39 elements composing the ground structure. The design vector v has thus 45 variables: v ¼ fh; LY ; LV ; x4 ; y4 ; x5 ; y5 g;

ð29Þ

where h is the vector with the thicknesses of the elements of the ground structure (he; e = 1, . . ., 39), LY and LV are the distance between the horizontal boundaries and length of the lateral borders of the parallelogram, and x4, y4 and x5, y5 are the position coordi-

s:t:

ðbÞ

qðvÞ ¼ qr  0:5 6 x4 =LH 6 0:5

ðcÞ

 0:5 6 y4 =LY 6 0:5  0:5 6 x5 =LH 6 0:5

ðdÞ ðeÞ

 0:5 6 y5 =LY 6 0:5

ðfÞ

LY =LV 6 1

ð32Þ

ðgÞ

and 0
ðhÞ;

where rle is the set of potentially critical stresses in all elements e of the ground structure for all loading conditions l, and where the relative density is given by T

qðvÞ ¼ h l=A

ð33Þ

with l containing the element lengths. Constraint (32b) fixes the relative density of the cellular material to a constant, qr. The internal nodes are constrained to remain inside the boundaries of the RVE by (32c–f). Relation (32g), limits the ratio LY/LV for obvious geometrical reasons. The thicknesses h are bounded in (32h) by a small positive number, k, in order to avoid singularity problems resulting from dangling nodes when the thicknesses of all the elements connected to the nodes approach zero. With regard to the objective function, we consider two types of failure modes: either the maximum axial stress in the element, rjAe , reaches the failure stress of the solid phase, rfs, assumed equal in

-x5 x4 8

H

+

10

9

V+

5

y5 LV

7

4

3

-y4

LY

n θn

Vy x

1

H-

2

6

LH Fig. 6. The ground structure, containing 39 rigid-jointed elements and 5 generic nodes (black dots), utilized for the strength design problem.

4022

F. Lipperman et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4016–4026

tension and compression, or, if the element has a compressive normal force, the maximum average compressive stress in the element, rCe, is equal to the elastic buckling critical stress, rbe, where rjAe

¼ rfs ;

ð34Þ

rCe ¼ rbe with rjAe ¼

jN e j 6jM je j þ ; 2 he he

rCe ¼  rbe ¼

j ¼ a; b;

ð35Þ

Ne ; he

ð36Þ

E 2 p ðhe =Le Þ2 : 12

ð37Þ

Here jNej is the absolute internal axial force in element e, and jMje j is the absolute value of the internal bending moment at extremity j, and Le is the length of element e. For conservative reasons the critical Euler buckling stress is calculated considering pin-jointed elements (see, for instance [27]). One will note that rbe decreases as the thickness of the element is reduced. This is very significant during the solution of a topology design problem, since including the buckling failure mode prevents the formation of slender elements and thus impedes the disappearance of non-effective elements. Therefore, following [28] the elastic buckling critical stress was modified rbe ¼

E 2 p ðhe =Le Þ2 expðh0 =he Þ; 12

ð38Þ

where h0 is a prescribed small thickness. To decide which failure mode is active one needs to relate E to rfs. We have assumed E/rfs = 103, which is a common ratio for ceramics and many metal materials [29]. Accordingly, the set of critical stresses is rle ¼

 l rCe raAe ; rbAe ; ce

ð39Þ

3 2

with ce ¼ 1012p ðhe =Le Þ2 expðh0 =he Þ. The design was carried out in two steps. In a first stage, (32) was solved using the h0 artifice (38). This is a standard topological procedure, augmented by some geometrical variables, where elements can reach the lower bound. In the next step the lower bound elements as well as all unconnected nodes were removed and problem (32) was solved again with topological variables only and including the nominal buckling (37). The actual minimization was performed with the fminimax subroutine of MATLAB which is a sequential quadratic programming algorithm.

σ'UA = {1, 0, 0}T

y'

y

5. Numerical results In this section we design two sets of two materials using the mixed topology/shape design formulation presented in Section 4. In a first instance two materials were designed: one for a uniaxial and one for a pure shear macroscopic state of stress, both of arbitrary orientations. In the next example one material was designed for a transversely isotropic tensional macroscopic stress field and the other for a transversely isotropic compressive macroscopic stress field. All cases were solved for 6 different relative densities: qr = 0.05, 0.1, 0.15, 0.2, 0.25, 0.3. 5.1. Uniaxial tension and pure shear of arbitrary orientation This paragraph considers two types of macroscopic stress fields, uniaxial and pure shear, applied in arbitrary orientations. In both cases, the orientation of the material is fixed in an x–y system of axes, while an x0 –y0 system, which is rotated by angle a, measured in a counterclockwise direction from axes x to x0 , defines the orientation of the applied macroscopic stresses (see Fig. 7). The macroscopic stress field in the x0 –y0 system of axes is denoted by r0 ¼ frx0 ; ry0 ; sx0 y0 gT . The applied uniaxial macroscopic stress field is denoted by r0UA ¼ f1; 0; 0gT and the pure shear macroscopic stress field by r0PS ¼ f0; 0; 1gT . To emulate the arbitrariness of the loading direction the problem was solved with multiple loading conditions where 50 sets of macroscopic stress fields were considered with orientations in the range 0 6 a 6 (49/50)p, with increments of p/50. In practice, the analysis of the cellular material in all directions is based on three analyses under unit macroscopic stress fields in the x–y coordinates. The design space had local optima and each optimization problem was therefore solved many times starting from different random initial designs. Interestingly, both the uniaxial and the pure shear macroscopic stress fields and for all considered relative densities converged consistently on one of three different topologies of triangulated patterns with almost identical strengths. These are the mixed square/triangular cells lattice, which is herein denoted by mc, the equilateral triangles (et) and the Kagomé lattice (ka) respectively shown in Fig. 8a–c. The optimal strengths, SUA and SPS, for the uniaxial and pure shear cases, as a functions of the relative density are given in Fig. 9. A perusal of the results will show that the mc and et grids fare in close vicinity (with a negligible advantage for mc) while the ka pattern exhibits a somewhat lower strength in the upper density range. This trend changes when low relative densities are considered where the Kagomé becomes the preferred lattice. The

σ'PS = {0, 0, 1}T

y' x' α

a

x

y x' α

x

b

Fig. 7. The orientation of the x0 –y0 system of axes for the (a) uniaxial and (b) pure shear macroscopic stress fields.

4023

F. Lipperman et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4016–4026

a

b

c

Fig. 8. The strength-optimal topologies for arbitrary uniaxial and pure shear macroscopic stress fields: (a) mixed cells lattice; (b) equilateral triangles lattice; (c) Kagomé lattice.

Table 3 The thickness to length ratios of the elements of the unit cells for arbitrary uniaxial or pure shear (in parenthesis when different from the uniaxial case) macroscopic stresses qr

0.05 0.10 0.15 0.20 0.25 0.30

(he/Le)mc Short elements

Long elements

0.020 0.050 (0.041) 0.075 0.100 0.125 0.150

0.015 0.025 (0.029) 0.037 0.050 0.063 0.075

(he/Le)et

(he/Le)ka

0.014 0.029 0.043 0.058 0.072 0.087

0.029 0.058 0.087 0.115 0.144 0.173

The values of the geometry design variables for the three optimal topologies are given in Table 2 (missing coordinates indicate removed nodes) while Table 3 details the thicknesses of the elements (topological variables). 5.2. Transversely isotropic macroscopic tension and compression

Fig. 9. The strength variation of the mixed cells, equilateral triangles and Kagomé lattices under (a) uniaxial and (b) pure shear macroscopic stresses, as a function of the relative density.

thickness to length ratio of the elements decreases with the relative density and, in consequence, the relative importance of the buckling failure mode increases. At such low relative densities, in the mc and et lattices, buckling overcomes the tensile rupture failure mode but not so for the Kagomé. (See [20], for a discussion on the benefits of the Kagomé grid).

Table 2 The geometry variables for arbitrary uniaxial and pure shear macroscopic stresses Topology

LY/LH

LV/LH

x4/LH

y4/LY

x5/LH

y5/LY

Mixed cells Equilateral triangles Kagomé

1 0.866 0.866

1 1 1

0 – 0

0 – 0

– – –

– – –

In this paragraph, the topology design problem was solved separately for two transversely isotropic macroscopic stress fields (single loading): tensional rTI = {1, 1, 0}T, and compressive rCI = {1, 1, 0}T. For the tensional case the optimal strengths, STI, obtained for the range of densities considered are shown by the straight line in Fig. 10b. As for the corresponding optimal topologies we find that for any given relative density the design space is literally full of different optimal patterns of same strength. We generated subsets of such topologies by starting from different initial designs. Some layouts, obtained for different densities, are shown in rows (a) and (c) of Fig. 10. The topologies in Fig. 10c, except the double hexagon layout (H–H), were known to have optimal strengths for adequately chosen thicknesses of the uniform elements (see, for instance, [8,30]). Less obvious are the patterns in Fig. 10a and the H–H grid. These are a few, out of a plethora of anisotropic cells which have optimal strength. It was easy to establish that all these optimal topologies, anisotropic as a rule, expand isotropically under macroscopic transversely isotropic tensional stresses. Indeed, all the elements have uniform extensional strain, even for intricate topologies. Such uniform internal strains imply a uniform macroscopic expansion of magnitude equal to the internal strains in the elements. Moreover, the effective compliance relations D computed for randomly chosen topologies had the property

4024

F. Lipperman et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4016–4026

a

b

S TI σ fs

ρr

ρr > 0.04

c

ρr > 0.07

ρr > 0.05

ρr > 0.06

ρr > 0.10

ρr > 0.12

Fig. 10. (b) Maximum strength as a function of relative density and corresponding optimal layouts of cellular materials under transversely isotropic macroscopic stress fields: (a) and (c) tensional stress for entire range of densities; (c) compressive stress, each layout valid for the indicated range below the figure.

D11 ¼ D22 ;

D13 ¼ D23 :

ð40Þ

For instance, the effective compliance matrix of the last cellular material shown in Fig. 10a, which optimal geometrical parameters and thicknesses are given in Table 4, are 8 9 9 0 18 16:9 > > < ex > = 1 16:8 3:5 < rx > = B C ey ¼ @ 3:5 16:8 16:9 A ry : ð41Þ > > > > E : ; : ; cyx syx 16:9 16:9 157:7 This is enough for isotropic expansion. A macroscopic stress field rx = 1, ry = 1, sxy = 0 yields indeed ex = ey, cxy = 0. One can show that for uniform expansion under isotropic tension Table 4 The geometrical parameters of the fourth strength optimal cellular material shown in Fig. 10a under a transversely isotropic macroscopic tensional stress field LY/LH

0.731

LV/LH

1.137

Element 1

Element 2

Element 3

U1 (rad)

h1/L1

U2 (rad)

h2/L2

U3 (rad)

h3/L3

0

4.32  102

1.396

8.95  102

2.443

1.32  102

STI qr ¼ ; rfs 2

ð42Þ

which is the relation that was established empirically in Fig. 10b. As for the transversely isotropic compressive case, it was noted that the optimal designs have shed all the elements which, if present, would have buckled. And since we have assumed a same failure stress in compression and in tension the optimal strength variation with relative density given in Fig. 10b is also valid for the compressive case. There is however a difference between the tension and compression cases. With compression, in order to avoid buckling, the material has a tendency to concentrate in fewer but bulkier elements. The designs in Fig. 10a do not appear in a compression field. Patterns which are strength optimal in compression are shown in Fig. 10c. And there is more. These patterns belong to the optimal set above some density threshold, indicated below each grid, below which they buckle and are therefore suboptimal. Interestingly, the triangle is optimal for qr > 0.12 whereas the hexagon stays optimal until qr = 0.04. Finally, it is easy to show that all these optimal designs have optimal effective bulk modulus. Indeed, borrowing from isotropic

F. Lipperman et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4016–4026

materials, the effective bulk modulus for any general topology was defined as the ratio between the average stress, r, and the dilatation, D j¼

r ðrx þ ry Þ=2 ¼ : D ex þ ey

ð43Þ

The macroscopic strain field of an optimal material in an isotropic macroscopic tensional stress field {STI, STI, 0}T is ex ¼ ey ¼

rfs 2STI ¼ ; E Eqr

ð44Þ

exy ¼ 0; which when substituted with (42) in (43) yields for this class of optimal materials j qr ¼ : E 4

ð45Þ

This value of j coincides with the Hashin–Shtrikman (H–S) upper bound [31] on the bulk modulus for 2D isotropic two-phase material in the low density range. There are published examples of media with effective isotropic or square symmetry which produce isotropic strain under isotropic stress and have maximum effective bulk modulus for the given density [8]. Herein we have shown that there exist a multitude of non-isotropic cellular materials without any symmetry which have uniform strain under macroscopic isotropic loading and hence reach the H–S bounds for isotropic materials. Furthermore, new optimal materials can be assembled by superimposing different topologies of materials which are independently optimal, as mentioned in [8]. The H–H layout in Fig. 10c (superposition of two irregular hexagon patterns) seems to be a point in case. Although the number of topologies corresponding to this class of materials is virtually unlimited, not for all layouts it is possible to find a thickness distribution such that under a macroscopic isotropic stress field a uniform macroscopic strain field is obtained. A typical example is the T–T cellular material analyzed earlier in this work. The equilibrium equations cannot be fulfilled with axial forces in the members.

4025

was emulated by using multiple loading conditions covering two quadrants and the maximum axial stress anywhere in the RVE for all loadings was minimized for a constant material density. The design space had several local minima and the procedure was therefore repeated with many randomly chosen initial designs. The algorithms converged systematically on one of the three following patterns: squares with two diagonals, equilateral triangles and Kagomé grids. These triangulated unit cells have pronounced axial modes and almost equal minimum stress with a slight advantage for the square in most of the considered range of densities. The following examples, transversely isotropic macroscopic tensile and compressive stress cases, were especially interesting since the design space is literally awash with local minima. It turned out that a plethora of very different geometries and topologies reached the same minimum stress condition for any given relative density. In contrast to the tension case, where very slender elements appear in the optimal layout, with compression such elements are removed due to the Euler buckling constraint. In consequence optimal layouts in the compressive case have fewer but bulkier elements. The common characteristic of all these optimum designs, both in tension and in compression, is that the materials have uniform strain under isotropic stresses although the unit cells were as a rule anisotropic. It was a small step to show that these optimum materials reach the Hashin–Shtrikman upper bound on the bulk modulus for isotropic materials in the low density range.

Acknowledgments The authors thank O. Sigmund and C.C. Seepersad for valuable information, Y. Benveniste for helpful discussions, P. Suquet for having brought to their attention the work of the French group, and an anonymous reviewer for having suggested to add the transversely isotropic macroscopic compression case. The authors also gratefully acknowledge the support of ISF, the Israel Science Foundation, in providing partial funding for the present research under Grant 838/06.

6. Summary In the context of optimal design of cellular material for maximum strength this paper has applied a unified method for stress and strain localization for periodic material where the unit cell is composed of rigid-jointed beam elements. Stress and strain periodic boundary conditions are applied to the RVE which leads to a linear system of equilibrium equations for the stiffness matrix of which a slightly modified assembly method is proposed. When a macroscopic stress field is imposed the equations will produce in one analysis the distribution of the axial and shear forces and bending moments in the elements of the RVE without passing through homogenization (three analyses). This reduction in the analysis effort is especially important in an optimization context where the optimum is reached after a large number of candidate designs have been analyzed. Moreover, the method can accommodate unit cells of any arbitrary topology and geometry. For the design of the RVE a ground-structure methodology was used where the topological design variables were the thicknesses of the elements. The latter were augmented by geometry variables controlling the boundaries of the parallelogram-shaped RVE and the coordinates of two internal nodes. Euler buckling of the uniform elements was included with a special modification to allow the elimination of slender elements. The technique was first employed for the design of optimal cellular materials subjected to a uniaxial stress fields and to pure shear stresses both of arbitrary orientations. The random orientation of the macroscopic stresses

References [1] P. Suquet, Elements of homogenization theory for inelastic solid mechanics, in: E. Sanchez-Palencia, A. Zaoui (Eds.), Homogenization Techniques for Composite Media, Lectures Notes in Physics, vol. 272, Springer-Verlag, 1987, pp. 193–278. [2] O. Debordes, C. Licht, J.J. Marigo, P. Mialon, J.C. Michel, P. Suquet, Limit loads of highly heterogeneous structures, in: J.P. Grellier, G.M. Campel (Eds.), Tendances actuelles en calcul des structures Pluralis, Paris, 1985, pp. 55–71 (in French). [3] J.C. Michel, H. Moulinec, P. Suquet, Effective properties of composite materials with periodic microstructure: a computational approach, Comput. Methods Appl. Mech. Engrg. 172 (1999) 109–143. [4] A. Anthoine, Derivation of the in-plane elastic characteristics of masonry through homogenization theory, Int. J. Solids Struct. 32 (1995) 137–163. [5] L.J. Gibson, M.F. Ashby, Cellular Solids: Structure and Properties, second ed., Cambridge Solid State Science Series, Cambridge University Press, Cambridge, UK, 1997. [6] O. Sigmund, Materials with prescribed constitutive parameters: an inverse homogenization problem, Int. J. Solids Struct. 31 (1994) 2313–2329. [7] O. Sigmund, Tailoring materials with prescribed elastic properties, Mech. Mater. 20 (1995) 351–368. [8] Z. Dimitrovova, L. Faria, New methodology to establish bounds on effective properties of cellular solids, Mech. Compos. Mater. Struct. 6 (1999) 331–346. [9] C. Hohe, W. Becker, Effective elastic properties of triangular grid structures, Compos. Struct. 45 (1999) 131–145. [10] J. Hohe, C. Beschorner, W. Becker, Effective elastic properties of hexagonal and quadrilateral grid structures, Compos. Struct. 46 (1999) 73–89. [11] Z. Dimitrovova, A new methodology to establish upper bounds on open-cell foam homogenized moduli, Struct. Multidisc. Optim. 29 (2005) 257–271. [12] C.C. Seepersad, J.K. Allen, D.L. McDowell, F. Mistree, Robust design of cellular materials with topological and dimensional imperfections, J. Mech. Des. 128 (2006) 1285–1297.

4026

F. Lipperman et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4016–4026

[13] J. Yan, G. Cheng, S. Liu, L. Liu, Comparison of prediction on effective elastic property and shape optimization of truss material with periodic microstructure, Int. J. Mech. Sci. 48 (2006) 400–413. [14] R.G. Hutchinson, N.A. Fleck, The structural performance of the periodic truss, J. Mech. Phys. Solids 54 (2006) 756–782. [15] L.J. Gibson, M.F. Ashby, G.S. Schajer, C.I. Robertson, The mechanics of twodimensional cellular materials, Proc. Roy. Soc. Lond. A 382 (1982) 25–42. [16] C. Chen, T.J. Lu, N.A. Fleck, Effect of imperfections on the yielding of twodimensional foams, J. Mech. Phys. Solids 47 (1999) 2235–2272. [17] A.-J. Wang, D.L. McDowell, In-plane stiffness and yield strength of periodic metal honeycombs, J. Engrg. Mater. Technol. 126 (2004) 137–156. [18] W.E. Warren, A.M. Kraynik, Foam mechanics: the linear elastic response of two-dimensional spatially periodic cellular materials, Mech. Mater. 6 (1987) 27–37. [19] D. Lukkassen, A. Meidell, S. Vigdergauz, On the elastic deformation of symmetric periodic structures, Quart. J. Mech. Appl. Math. 56 (2003) 441–454. [20] S. Hyun, S. Torquato, Optimal and manufacturable two-dimensional, Kagomélike cellular materials, J. Mater. Res. 17 (2002) 137–144. [21] J.M. Guedes, H.C. Rodrigues, M.P. Bendsoe, A material optimization model to approximate energy bounds for cellular materials under multiload conditions, Struct. Multidisc. Optim. 25 (2003) 446–452.

[22] W.S. Dorn, R.E. Gomory, H.J. Greenberg, Automatic design of optimal structures, J. Méchanique 3 (1964) 25–52. [23] M.W. Dobbs, L.P. Felton, Optimization of truss geometry, J. Struct. Div. ASCE 95 ST10 (1969) 2105–2118. [24] M. Zhou, G.I.N. Rozvany, DCOC – An optimality criteria method for large systems. 1. Theory, Struct. Optim. 5 (1992) 12–25. [25] S.D. Guest, J.W. Hutchinson, On the determinacy of repetitive structures, J. Mech. Phys. Solids 51 (2003) 383–391. [26] J. Renton, On the analysis of triangular mesh grillages, Int. J. Solids Struct. 2 (1966) 307–318. [27] A.M. Hayes, A. Wang, B.M. Dempsey, D.L. McDowell, Mechanics of linear cellular alloys, Mech. Mater. 36 (2004) 691–713. [28] G.I.N. Rozvany, Difficulties in truss topology optimization with stress, local buckling and system stability constraints, Struct. Optim. 11 (1996) 213–217. [29] W. Bolton, Engineering Materials Technology, third ed., ButterworthHeinemann, Oxford, 1998. [30] S. Torquato, L.V. Gibiansky, M.J. Silva, L.J. Gibson, Effective mechanical and transport properties of cellular solids, Int. J. Mech. Sci. 40 (1998) 71–82. [31] Z. Hashin, S. Shtrikman, A variational approach to the theory of the elastic behaviour of multiphase materials, J. Mech. Phys. Solids 11 (1963) 127–140.