Sandwich plates of minimal compliance

Comput. Methods Appl. Mech. Engrg. 197 (2008) 4866–4881

Contents lists available at ScienceDirect

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

Sandwich plates of minimal compliance S. Czarnecki a, M. Kursa b, T. Lewin´ski a,* a b

Institute of Structural Mechanics, Warsaw University of Technology, Al. Armii Ludowej 16, 00-637 Warsaw, Poland Institute of Fundamental Technological Research, Polish Academy of Sciences, Ul. S´wießtokrzyska 21, 00-049 Warsaw, Poland

a r t i c l e

i n f o

Article history: Received 11 April 2008 Received in revised form 30 June 2008 Accepted 17 July 2008 Available online 5 August 2008 Keywords: Minimum compliance problem Sandwich plates Topology optimization

a b s t r a c t The subject of the paper is an optimal choice of material parameters characterizing the core layer of sandwich plates within the framework of the conventional plate theory in which the core layer is treated as soft in the in-plane direction. The mathematical description is similar to the Hencky–Reissner model of plates with transverse shear deformation. Here, however, the bending stiffnesses and the transverse shear stiffnesses can be designed independently. The present paper deals only with optimal design of the core layer to make the plate compliance minimal. Two core materials are at our disposal, which leads to the ill-posed problem. To consider it one should relax this problem by admitting composite domains and characterize their overall properties by the homogenization formulae. The numerical approach is based on this relaxed formulation thus making it mesh-independent. The equilibrium problem is solved by the DSG3 finite element method. The optimization results are found with using the convergent updating schemes of the COC method. Crown Copyright Ó 2008 Published by Elsevier B.V. All rights reserved.

1. Introduction The minimum compliance problem for the bodies composed of two constituents of isotropic and linearly elastic characteristics, with the isoperimetric condition imposed on the amount of one constituent, was solved in 2D case by Gibiansky and Cherkaev [1], see Cherkaev [2] and Allaire [3]. The initial problem of forming the body from two constituents is ill-posed in both 2D and 3D settings and should be replaced by its relaxation in which the body is endowed with a composite microstructure, see Lipton [4]. The main theorem states that the minimum of the compliance is achieved by so called hierarchical stiff laminates of third-rank in the 3D case and of second-rank in the 2D case. If one assumes that the two materials at our disposal are ordered (this means that j2 > j1 and l2 > l1 , where ja , la are the bulk and shear moduli, a ¼ 1; 2), then the stiff laminates of third- or second-rank should be constructed by a subsequent layering of the stronger material with the previously obtained laminate, starting point being a conventional, first-rank two-phase laminate. Thus, the softer material is used only at the first step of the layering process. Although many authors used this result in 3D analyses in 1980s and 1990s, the only explicit proof of this theorem in 3D case known to the present authors is that reported by Allaire [3]. The numerical procedures for the 3D case were developed by Jacobsen et al. [5], Borrvall ´ az and Lipton [7], Olhoff et al. [8] and Czarand Petersson [6], Dı necki and Lewin´ski [9]. The relaxed formulation makes it possible

* Corresponding author. Tel.: +48 22 660 6379; fax: +48 22 825 69 85. E-mail address: [email protected] (T. Lewin´ski).

to consider the case of the softer material being the void, which means the generalized shape design. In the 2D case this relaxed formulation is degenerated while is well posed in the 3D case, see Allaire and Kohn [10] and Allaire et al. [11]. The numerical methods of solving the 2D and 3D relaxed problems and their approximations are discussed in Bendsøe [12] and Bendsøe and Sigmund [13,14]. The relevant minimum compliance problem for the thin plates in bending was solved a little earlier than 2D elasticity problem, namely in Gibiansky and Cherkaev [15]. The optimization problem can be reduced to an equilibrium problem of a hyperelastic plate, see Lewin´ski and Telega [16,Chapter 26]. The numerical methods are discussed in Bendsøe [12], Krog and Olhoff [17], Czarnecki _ and Lewin´ski [18], Kolanek and Lewin´ski [19] and Dzierzanowski [20]. We do not mention here the papers on optimization of thin plates which were based on the initial and not relaxed optimal layout formulation. The minimum compliance problem of plates loaded simultaneously in-plane and out-of-plane is similar to the optimum design of shells, especially-shallow shells. These problems involve two stress fields, hence are coupled, which introduces new essential difficulties. The relaxed formulation of these problems are dis_ cussed in Dzierzanowski [20] and Lewin´ski [21]. The optimization problem is there reduced to an equilibrium problem of a hyperelastic shell or plate. The minimum compliance problem of two component plates of moderate thickness with transversely distributed material proper´ az et al. [22], where the so ties was the subject of the paper by Dı called moment formulation was used, see Francfort et al. [23], to

0045-7825/$ - see front matter Crown Copyright Ó 2008 Published by Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2008.07.005

S. Czarnecki et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4866–4881

model the underlying third-rank microstructure of the effective ´ az et al. [22] is probably the composite plate. Just this paper by Dı closest in spirit to the present paper amongst all the hitherto existing papers on the optimal plate design: it is based on the same plate theory with transverse shear deformation (yet applied to a transversely homogeneous plate), it is also aimed at minimizing the compliance of the two-component plate and its main subject is also focused on the appropriate reformulation of the initial two-component problem into the composite-like problem with an underlying microstructure. The third-rank composites used by ´ az et al. [22] are sufficient to realize the minimum of the compliDı ance. In this manner the authors have circumvented the still open problem of full relaxation of the minimum compliance design of moderately thick two component plates. The bending and transverse shearing deformations of sandwich plates of soft core can be described by the conventional theory of plates of moderate thickness, as noted by Reissner [24]. The modeling is conditioned by possibility of neglecting the own bending stiffnesses of the faces and reducing the core layer deformations to the transverse direction. This theory can be viewed as a simplification of the Hoff’s theory, see Lewin´ski [25]. The faces of the sandwich plates are usually of constant thickness. In the present paper, their characteristics will not be designed. The distance between the faces will also be treated as fixed. Similarly as in the papers by Buannic et al. [26] the core layer is subject to the designing. This layer keeps distance between the faces, being very soft inplane. One of possible design is the honeycomb core. It can be treated as a transversely isotropic or transversely orthotropic material, depending on the regularity of its in-plane geometry. Within such a modelling the effective transverse moduli G13 , G23 of the core layer depend on the elastic moduli and dimensions of the honeycomb walls, the appropriate formulae being given in Gibson and Ashby [27] and Shi and Tong [28]. The transverse stiffnesses Hab of the sandwich plate are expressed by conventional formulae involving the moduli G13 , G23 as well as transverse dimensions of the faces and the core. This method of computing the stiffnesses Hab will be assumed in the present paper, although this method neglects that these stiffnesses depend on the thickness of the core in a much more complex manner, as was indicated by Becker [29]. The effective stiffnesses of sandwich plates should be assessed by the methods which treat the faces and the core as one body. Saying this in the terms of homogenization theory, the basic cell of a plate of repetitive structure is a three-dimensional body and its 3D analysis should determine the effective stiffnesses. Such approach is proposed by Hohe [30], Buannic et al. [26]. However, the 3D approach makes the whole optimization process very complex; it will not be applied here to make the optimization as simple as possible. The subject of the paper will be an optimal layout of two materials of the core of the honeycomb structure to minimize the compliance of the sandwich plate of fixed faces, of fixed distance between the faces and of fixed in-plane dimensions. The optimization makes the plate as stiff as possible with respect to the transverse loading. The layout of the characteristics of the core layer

4867

adapts itself to the distribution of the transverse shear forces and to the transverse shear deformations. The relaxation technique applied is based on the results of the homogenization formulae for the effective stiffnesses of the sandwich plate of periodic distribution of transverse shear characteristics. The so called in-plane scaling approach is used, see Bourgeat and Tapiero [31], Lewin´ski [32], Lewin´ski and Telega [16,Chapter 5.2]. The basic cell problem decouples into two problems: for bending and transverse shearing. The former problem has trivial solution, while the mathematical structure of the latter is relatively simple, since it is a scalar elliptic problem with periodic coefficients, known from the homogenization theory of the temperature distribution problem. This analogy makes it possible to transfer the known theoretical homogenization results and optimal layout results from the scalar elliptic problem (see Allaire [3], Chapters 1 and 3) to the problem of optimal design of a core of a sandwich plate. In the COC method to be used, see Bendsøe and Sigmund [13], the linear equilibrium problem must be solved many times, hence a reliable FEM approach for sandwich plate is necessary. We shall use here the DSG3 triangular element developed by Bletzinger et al. [33] for Hencky plates. This element passes the consistency test of Rakowski [34], as proved in Kursa and Lewin´ski [35]. The test means that in the case of regular mesh the FEM equations assume the form of the recurrent algebraic equations which converge (when the size of the triangles tend to zero) to the three partial differential equations describing the local equilibrium equations of the plate. However, this element is sensitive to the nodes numeration, which is very annoying. The numeration has an influence on the convergence speed but in all the cases of numeration considered the convergence held good. The following conventions are adopted: Latin indices run over 1, 2, 3, while the Greek indices assume the values 1, 2. The summation convention refers to the repeated indices at different levels. 2. Eric Reissner model of sandwich plate of soft core Consider a plate composed of the three layers of thicknesses d, 2c, d, symmetric with respect to the middle plane X of the core layer. The left-handed Cartesian system ðx1 ; x2 ; x3 Þ is adopted, while ðx1 ; x2 Þ parameterize the domain X and jx3 j < c þ d (see Fig. 1). The faces and the core are glued strongly such that the displacements ðw1 ; w2 ; w3 Þ within the plate can be viewed as continuous functions; wi is directed along the xi axis. The material of the faces is elastic with moduli Cf ¼ ðC ijkl f Þ; the x3 axis is assumed as the axis a ¼ 0, C 3f abc ¼ 0, a; b; c ¼ 1; 2. Let of material symmetry or C 333 f

e abkl ¼ C abkl  C ab33 C 33kl =C 3333 C f f f f f

ð2:1Þ

be components of the tensor of reduced moduli of the generalized plane stress state. The bending stiffnesses are assumed according to the formula

Fig. 1. Sandwich plate.

4868

S. Czarnecki et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4866–4881

2 e abkl Dabkl ¼ 2db C ; f

ð2:2Þ

where b ¼ c þ d=2. This formula can be derived by imposing the Hoff type kinematics (see Fig. 6.1.2 in Lewin´ski and Telega [16]) and then neglecting the own bending stiffnesses of the faces in the Hoff plate model, see Sections 6.1 and 6.3 in Ref. [16]. The transverse shear stiffnesses are in this model determined by the transverse shear moduli of the core. If the Hooke tensor of the then the stiffnesses due to core material has components C ijkl c transverse shear depend only on C ca3b3 . We assume that

C 1313 ¼ C 2323 ; c c

C 1323 ¼0 c

ð2:3Þ

can be considered independently of the in-plane deformation. The deformation of the plate is determined by the angles of rotation of the plate cross-sections u1 ðx1 ; x2 Þ, u2 ðx1 ; x2 Þ and the deflection w3 ¼ wðx1 ; x2 Þ, see the sign conventions in Fig. 3. The deformations are defined as in the Hencky–Reissner plate theory

qab ¼ eab ðuÞ; ca ¼ ba ðu; wÞ;

where u ¼ ðu1 ; u2 Þ and the operators eab ðÞ and ba ð; Þ are defined as below

eab ðuÞ ¼

and write

C a3b3 ¼ lc dab ;

ð2:4Þ

where lc is the modulus of shear. The stiffnesses due to transverse shear of the plate are given by

Hab ¼ Hdab ;

H¼k

2b c

2

lc ;

ð2:5Þ

which can be derived by a consistent derivation of the Hoff theory, see Section 6.1 in [16]. The shear correction factor k will be assumed as k ¼ 1. The same formulae (2.2) and (2.5) hold in the simplified Reissner theory [24] which will be used in this work. The moduli e abkl and l are chosen independently. In the theory considered C c f the elastic moduli of the core do not affect the bending stiffnesses abkl while the elastic moduli of the faces C ijkl do not affect the D f transverse shear stiffnesses Hab of the plate. The core material is usually cellular and its main role is to keep the distance 2c constant. This material can be weak in the in-plane directions x1 and x2 . In particular, a honeycomb core can be applied, see Fig. 2. If h ¼ p=6 and h ¼ l, t 1 ¼ t 2 ¼ t the core has isotropic properties ¼ C 2323 ¼ lc can be assessed by the formula and the moduli C 1313 c c

lc

pffiffiffi 3t ¼ l; 3 l

ð2:6Þ

where l is the shear modulus of the material of the walls of the core, see Gibson and Ashby [27]. The honeycomb core is usually produced by stretching the metal sheets, hence some walls have double thickness, e.g. t 1 ¼ 2t, –C 2323 . This effect will not be taken into t2 ¼ t. Consequently, C 1313 c c account to keep the isotropy of the formula (2.5). Eqs. (2.5) and (2.6) are strongly simplified. More accurate formulae for the moduli Hab can be found in Shi and Tong [28] and in Hohe [30], Becker [29]. We shall assume that the plate is subject to a transverse loading of intensity pðx1 ; x2 Þ such that the bending and transverse shearing

ð2:7Þ

1 ðu þ ub;a Þ; 2 a;b

ba ðu; wÞ ¼ ua þ w;a

ð2:8Þ

and ðÞ;a ¼ o=oxa . We define the stress resultants, the moments and the shear forces

Mab ¼ Dabkl qkl : Q a ¼ Hab cb :

ð2:9Þ

Let ðn; sÞ be unit vectors outward normal and tangent to the boundary line C ¼ oX

n ¼ ðn1 ; n2 Þ;

s ¼ ðn2 ; n1 Þ;

n1 ¼ cos a;

n2 ¼ sin a;

a ¼ ðe1 ; nÞ; where ðe1 ; e2 Þ is the orthogonal basis of the system ðx1 ; x2 Þ. The boundary conditions concern

M n or un ;

M s or us ;

ð2:10Þ

Q n or w; where

Mn ¼ M ab na nb ;

M s ¼ M ab nb sa

Q n ¼ Q a na are the stress resultants referred to the basis ðn; sÞ and un ¼ u  n, us ¼ u  s are the angles of rotation of directions tangent and normal to the boundary C. In particular, along the simple support we have: w ¼ 0, M n ¼ 0, M s ¼ 0, while along the hinge support there holds: w ¼ 0, M n ¼ 0, us ¼ 0. Let U be the set of admissible displacements u ¼ ðu1 ; u2 ; wÞ satisfying the homogeneous kinematic boundary conditions, e.g. conditions of clamping on C0 , the part of the boundary C. We shall not consider the non-homogeneous kinematic boundary conditions to make use of the set U being the linear space. The equilibrium conditions of the plate have the form of the variational equation

Z

½M ab eab ðwÞ þ Q a ba ðw; vÞdx ¼ f ðw; vÞ;

ð2:11Þ

X

where w ¼ ðw1 ; w2 Þ are the trial angles of rotation and v represents the trial deflection. The right-hand side of (2.11) or f ðw; vÞ repre-

Fig. 2. Honeycomb core.

Fig. 3. Sign conventions.

4869

S. Czarnecki et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4866–4881

sents the virtual work of the loading. In the case of the transverse loading pðx1 ; x2 Þ this linear form is

f ðvÞ ¼

Z

pvdx:

ð2:12Þ

X

Substitution of (2.9) into (2.11) leads to the variational equation having unique solution, due to symmetry, ellipticity and continuity of the bilinear form which appears at the left hand side, putting aside subtle regularity problems of the functions from the set U, cf. Section 7.2 in Lewin´ski [25]. The trial stress resultant fields of moments N ¼ ðN ab Þ and transverse shear T ¼ ðT a Þ are statically admissible if they satisfy the variational Eq. (2.11). The set of such fields R is linear affine. Let us define the functionals

Jðw; vÞ ¼

1 2

Z

½Dabkl eab ðwÞekl ðwÞ þ Hak ba ðw; vÞbk ðw; vÞdx  f ðvÞ

X

5.2 in Lewin´ski and Telega [16], referred to the so-called in-plane scaling. The first homogenization results of moderately thick plates, based on this scaling, were found by Bourgeat and Tapiero [31].  The effective compliances ðhab Þ constitute the tensor inverse to   H , or h ¼ H1 . Tensor h can be determined directly from the formula

Q a hab Q b ¼ minfhq : ðhðyÞqÞijq 2 Rper ðYÞ; hqi ¼ Q g; 

where hi means averaging over Y and

Rper ðYÞ ¼ fq ¼ ðq1 ; q2 Þjdivq ¼ 0 and q  m assume opposite values at opposite sides of Yg

( ab

CðN; TÞ ¼

ab

H ðy1 ; y2 Þ ¼ H ðy2 Þ ¼

Z

ab

ðdabkl N N

kl

a b

þ hab T T Þdx;

ð2:14Þ

ð3:5Þ

and qa are appropriately regular, while m is the unit vector outward normal to oY. In the sequel, the strip-wise plates will be important, of the basic cell Y divided into two strips, as in Fig. 4. The effective transverse shear stiffnesses of such plate are given by

ð2:13Þ and

ð3:4Þ

H2 dab

for y2 2 ð0; hl2 Þ;

ab

for y2 2 ðhl2 ; l2 Þ;

H1 d

ð3:6Þ

where d ¼ D1 , h ¼ H1 ; these tensors appear in the equations inverse to (2.9)

where h 2 ð0; 1Þ; H1 , H2 are given positive constants. The formulae for the effective stiffnesses can be determined by Tartar’s rule, see [16, Eq. 5.6.42]

qab ¼ dabkl Mkl ; cb ¼ hba Q a :

h1 ðH2  H Þ1 ¼ ðH2  H1 Þ1 

X

ð2:15Þ

One can prove that the solution ðu1 ; u2 ; wÞ is minimizer of the functional J

Jðu; wÞ ¼ minfJðw; vÞjðw; vÞ 2 Ug

ð2:16Þ

and the pair ðM; Q Þ is minimizer of the functional (2.14)

CðM; Q Þ ¼ minfCðN; TÞjðN; TÞ 2 Rg:

h2 na H2ab nb

H2 ¼ H2 I;

I ¼ ðdab Þ:

H1 ¼ H1 I;

ð2:17Þ

On the basis of (3.7) and (3.8) we find the components of H

2

h1 H1 þ h2 H2



0

3 1 5:

3. Homogenization of transverse shear stiffnesses

h ¼

Let us consider a family of problems (2.16) indexed by a parameter e, hence now J ¼ J e and

where

H ¼4 

X

½Dabkl eab ðwÞekl ðwÞ þ Heak ðxÞba ðw; vÞbk ðw; vÞdx  f ðvÞ; ð3:1Þ

ab

ab

ab

where He ðxÞ ¼ H ðx=eÞ; x ¼ ðx1 ; x2 Þ, the quantities H ðÞ are Yperiodic, Y being a given rectangle ð0; l1 Þ  ð0; l2 Þ. The family of problems defined this way is denoted by ðPe Þ. The stiffnesses Hae b ðxÞ are eY-periodic. Let ðue ; we Þ be the solution of the problem ðPe Þ. One can prove that this solution converges weakly in a Sobolev space to the solution ðu ; w Þ of the homogenized problem ðP Þ in the form (2.16) with the potential in which the tensor H is replaced with its homogenized counterpart H . The effective tensor H is determined by the formula:

na Hab nb ¼ min

u2H1per ðYÞ

Z Y

!   ou ou dy; Hab ðyÞ na þ nb þ oya oyb

h1 H1

0

þ Hh22



hðhÞ 0

0



;  hðhÞ

hðhÞ ¼ ½ð1  hÞH1 þ hH2 1 ;  hðhÞ ¼ ð1  hÞh1 þ hh2 :

ð3:10Þ

ð3:11Þ

The plate of the layout of Fig. 4 has extremal effective properties, which will be used in the minimum compliance problem. 4. The minimum compliance problem of a sandwich plate with a composite core Assume that two materials to fill up the core layer are at our disposal. They differ in the modulus lc . Thus the transverse shear stiffnesses are given as H1ab and H2ab where

ð3:2Þ

where y ¼ ðy1 ; y2 Þ 2 Y and

H1per ðYÞ ¼ fv 2 H1 ðYÞjv assumes equal values at opposite sides of Yg: ð3:3Þ The formula (3.2) is the same as in the theory of homogenization of the conductivity problem, see Eq. (1.14) in Allaire [3]. The result (3.2) follows from the assumption De ¼ D in the analysis of Section

ð3:9Þ





1 2

ð3:8Þ 

Let h ¼ H1 , ha ¼ 1=Ha . Hence

J e ðw; vÞ ¼

ð3:7Þ

where h1 ¼ 1  h2 , h2 ¼ h and n ¼ ðn1 ; n2 Þ is a unit vector orthogonal to the strips; here n ¼ e2 ¼ ð0; 1Þ. Moreover

We omit here the regularity issues which decide whether the minimizers exist. The properties (2.16) and (2.17) can be proved by the methods of Necˇas and Hlavacˇek [36], applied there to the equilibrium problem of linear elasticity.

Z

n  n;

Fig. 4. Strip-like microstructure. Basic cell Y.

4870

Harb ¼

S. Czarnecki et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4866–4881

2b c

2

lcðrÞ dab ; r ¼ 1; 2:

ð4:1Þ

We shall assume that the core materials are of honeycomb structure, made from two different materials of various moduli l1 and l2 and equal dimensions t and l or

pffiffiffi 3t lðcrÞ ¼ l : 3 l r

ð4:2Þ

We shall assume that

Hr ¼

2b c

pffiffiffi 3t l ; 3 l r

Hr ¼ H r I

Z ðPk Þ : min CðvÞ þ k

ð4:3Þ

X

b k Þ : min ðP

ð4:4Þ

and the distribution of the stiffness values is given the same way

HðxÞ ¼ H1 v1 ðxÞ þ H2 v2 ðxÞ:

ð4:5Þ

SðTÞ ¼

ð4:7Þ

where

Cðv2 ; N; TÞ ¼

X

½dabkl Nab Nkl þ hab ðv2 ÞT a T b dx:

ð4:8Þ

Here, the matrix ðhab Þ is inverse to ðHab Þ, or h ¼ H1 and

hðxÞ ¼ h1 v1 ðxÞ þ h2 v2 ðxÞ; hr ¼ hr I;

hðxÞ ¼ h1 v1 ðxÞ þ h2 v2 ðxÞ

hr ¼ ðHr Þ1

ð4:9Þ

and we have h2 < h1 . The function v2 is bounded in X and assumes values 0 and 1 which is written as v2 2 L1 ðX; f0; 1gÞ. The amount of both the materials is viewed as given. Thus the volume of the cylindrical domains Xa  ðc  d; c þ dÞ, a ¼ 1; 2, is given. It is sufficient to impose the condition jX2 j ¼ A2 or

Z

X

v2 dx ¼ A2

ð4:10Þ

Z

v2L ðX;f0;1gÞ

(

W k ðT; vÞdx

ð4:12Þ

X

h1 kTk2

if

h2 kTk2 þ k if

v ¼ 0; v ¼ 1:

ð4:13Þ

Here, kTk2 ¼ ðT 1 Þ2 þ ðT 2 Þ2 . We write the functional SðTÞ in the form

or

Cðv2 Þ ¼ minfCðv2 ; N; TÞjðN; TÞ 2 Rg;

min 1

W k ðT; vÞ ¼

Cðv2 Þ ¼ f ðwÞ: According to (2.17) one can write

X

and

SðTÞ ¼

ð4:6Þ

 Nab dabkl Nkl dx þ SðTÞ ;

where

The function v2 is the design variable which determines the transverse stiffness of the plate and, consequently, all the fields ua , w, Mab , Q a describing its deformation and stress resultants. These relations will not be put explicit, we shall write only that C ¼ Cðv2 Þ, where C represents the plate compliance

Z

Z

ðN;TÞ2R

with r ¼ 1; 2, where I represents the unit tensor of second-rank. Note that H2 > H1 . Let X1 , X2 be subdomains of X around which both the core materials are placed, contributing to the stiffnesses H1 and H2 . Let va be a characteristic function of the domain Xa or va ðxÞ ¼ 1 if x 2 Xa and va ðxÞ ¼ 0 if x 2 X n Xa . The transverse stiffness tensor is of the form

HðxÞ ¼ H1 v1 ðxÞ þ H2 v2 ðxÞ



vdxjv 2 L1 ðX; f0; 1gÞ ;

where k is the Lagrangian multiplier. The problem ðPk Þ is re-written in the form

l2 > l1 . By (4.1) and (4.2) we have 2

Harb ¼ Hr dab ;

convexified and the vector problem-quasiconvexified. It occurs that the problem of optimal design of the transverse stiffness can be made well posed by convexification, as will be shown below. To show necessity of the relaxation let us replace the ðPÞ problem with

Z

U k ðT; xÞdx

ð4:14Þ

X

with

U k ðT; vÞ ¼ minfh1 kTk2 ; h2 kTk2 þ kg

ð4:15Þ

U k ðT; vÞ ¼ uk ðkTðxÞkÞ;

ð4:16Þ

uk ðzÞ ¼ minfh1 z2 ; h2 z2 þ kg:

ð4:17Þ

The function z ! uk ðzÞ is non-convex, see Fig. 5; we remember that h1 > h2 . The multiplier k > 0 if the condition (4.10) is bounded from below. The non-convexity of uk ðzÞ is the source of theoretical problems with ðP k Þ. These theoretical issues emerge as ill-posed numerical schemes suffering from mesh-dependency of the approximate solutions. 5. Relaxation of problem ðP k Þ Consider the minimizing sequences of layouts of both the materials, given by the sequences vn , such that the subsequent layouts vnþ1 contribute to the correction of the compliance, or Cðvnþ1 Þ < Cðvn Þ. This correction is usually achieved by making

to make the optimization problem reasonable. Of all possible layouts of both the core materials satisfying (4.10) we choose the layout realizing the minimum of the compliance. This problem is written as

min CðvÞ

Z



vdx ¼ A2 ; v 2 L1 ðX; f0; 1gÞ ;

ð4:11Þ

X

where v ¼ v2 to make further notation as simple as possible. The problem above will be named ðPÞ and it is the minimum compliance problem in its initial formulation. The problem ðPÞ has much in common with the layout problems in the conductivity theory involving the scalar temperature field and in the elasticity theory. The latter problems need relaxation, cf. Allaire [3, Sections 3.1, and 4.1], and a similar reformulation should be applied to the initial formulation ðPÞ. It occurs that the scalar problem should be

Fig. 5. Juxtaposition of two parabolas shows non-convexity of uk ðzÞ.

4871

S. Czarnecki et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4866–4881

the layout even more ground down, fragmented into smaller pieces. A minimizing sequence ðvn Þ determines the sequences of stiffness tensors Hn ðxÞ and compliance tensors hn ðxÞ which then determine the moments Mn and transverse forces Q n . From the sequence vn one can extract a subsequence weak-* converging to certain function m2 2 L1 ðX; ½0; 1Þ which represents the density of the stronger material No. 2. This function determines a generalized layout referring to a composite two component material of volume (or rather – area) fractions m1 ¼ 1  m2 and m2 . The mathematical aspects of the passage to the relaxed problem in the conductivity and elasticity settings can be found in Allaire [3]. The same arguments apply in the problem considered. Therefore, the mathematical subtleties will be omitted to concentrate on the very formulation of the relaxed problem. We shall draw upon the known results on homogenization of conductivity moduli of a periodic composite and the transverse stiffnesses of moderately thick plates of periodic structure. The clue is the mathematical analogy between the basic cell problems of homogenization of conductivity moduli and transverse stiffnesses of plates, cf. Section 1.1.3 in Allaire [3] and Section 5.2 in Lewin´ski and Telega [16]. Concluding, the optimal solution will refer to a plate with a non-homogeneous composite core with non-uniform effective transverse stiffnesses ðHab ðxÞÞ, x 2 X. The stiffnesses Hab ðxÞ are locally characterized by a periodic structure. Let us imagine an eYperiodic plate, Y being a rectangle parametrized by ðy1 ; y2 Þ in which the layouts of both the materials are determined by the function y ! vY2 ðx; yÞ, y ¼ ðy1 ; y2 Þ 2 Y. The distribution of transverse stiffnesses within the basic cell is

Hab ðx; yÞ ¼ Ha1b vY1 ðx; yÞ þ H2ab vY2 ðx; yÞ;

ð5:1Þ

where vY1 ðx; yÞ ¼ 1  vY2 ðx; yÞ. Similarly

hðx; yÞ ¼ h1 vY1 ðx; yÞ þ h2 vY2 ðx; yÞ:

Z Y

vY2 ðx; yÞdy ¼ m2 ðxÞ:

q 2 Rper ðYÞ; hqi ¼ T; hvY i ¼ h; vY 2 L1 ðY; f0; 1gÞg:

ð5:2Þ

ð5:3Þ

ð5:4Þ

1

hq : hðyÞqi P hqi : hh ðyÞi1 hqi;

ð5:9Þ

see Section 21.5 in Lewin´ski and Telega [16]. Hence, we estimate 1

hq : hðyÞqi P T : ðhh ðyÞi1 TÞ:

ð5:10Þ

1

According to (5.1) h ðyÞ ¼ HðyÞ or 1

hh ðyÞi ¼ hð1  vðyÞÞH1 þ vðyÞH2 i ¼ ð1  hÞH1 þ hH2 :

ð5:11Þ

According to (4.1) we assume that the tensors H1 and H2 are isotropic of the form (4.3). Thus, we write 1

hh ðyÞi1 ¼ ½ð1  hÞH1 þ hH2 1 I

ð5:12Þ

and the estimate (5.10) assumes the form

hq : ½ð1  vY ðyÞÞh1 þ vY ðyÞh2 qi P hðhÞkTk2 ;

ð5:13Þ

where hðhÞ is given by (3.11). If we indicate v realizing the equality in (5.13), the minimum problem (5.8) will be solved thus giving the explicit form of W  ðT; hÞ. The expression for W  is determined by the effective complianc  es hab , see (5.6). If one assumes vY ¼ vY2 as in Fig. 4, the tensor h will have the form (3.10). Let us direct the axes y1 and y2 of the cell Y (Fig. 4) along T such that T be directed along y1 . Then T 1 ¼ kTk, T 2 ¼ 0 and Y



ð5:14Þ

which proves that the above choice of the direction of ðy1 ; y2 Þ with respect to T and the choice of the division of Y into Y 1 and Y 2 realizes equality in (5.13) and leads to minimum in (5.6). Thus

2W  ðT; hÞ ¼ hðhÞkTk2 :

ð5:15Þ

The optimal plate has a strip-wise microstructure with the strips directed along the maximal shear. In the direction orthogonal to the strips the shear is zero. The vector field ^c ¼ ðc2 ; c1 Þ shows the transverse shear better than c ¼ ðc1 ; c2 Þ, see Fig. 6. Let us come back to (5.5). We put the compliance in the form below

e 2 ; N; TÞ ¼ Cðm

Z

dabkl Nab Nkl dx þ

X

Z

hðm2 ÞkTk2 dx:

ð5:16Þ

X

e (5.7), assumes the form The problem ð PÞ,

e ðm2 Þ þ k minf C

and

Z

m2 ðxÞdx j m2 2 L1 ðX; ½0; 1Þg;

ð5:17Þ

X

e 2 ; N; TÞ ¼ Cðm

Z

½dabkl Nab Nkl þ 2W  ðT; m2 Þdx:

ð5:5Þ

X

Let h 2 ½0; 1. The integrand W  ðT; hÞ is defined as below

2W  ðT; hÞ ¼ minfT a hab T b j vY 2 L1 ðY; f0; 1gÞ; hvY i ¼ hg; 

with

ð5:8Þ

The minimized expression hq : hðyÞqi can be estimated from below according to



The left hand side of the above equation will be denoted by hv2 ðx; Þi where hi is the averaging operator over Y. If the sequence vn tends weak-* to m2 then the sequence Cðvn Þ given by (4.7), (4.8) tends to e ðm2 Þ, where C

e ðm2 Þ ¼ minf Cðm e 2 ; N; TÞjðN; TÞ 2 Rg C

2W  ðT; hÞ ¼ minfhq : ½ð1  vY ðyÞÞ=½h1 þ vY ðyÞÞh2 qij

T a hab T b ¼ ðT 1 Þ2 h11 ¼ hðhÞkTk2 ;

The distribution of stiffnesses (5.1) within Y determines the effective stiffnesses ðHab ðxÞÞ according to the formula (3.2); this distribution does not affect the bending stiffnesses Dabkl which are assumed as uniform in X. The area fraction m2 ðxÞ determines the amount of the material 2 within the cell Y, referred to a given point x, or

1 jYj

Let us put (3.4) and (5.6) together

min ð5:6Þ

vY ¼ vY2 . The problem (4.11) is replaced by

Z  e ðm2 Þ m2 dx ¼ A2 ; m2 2 L1 ½X; ½0; 1 : e : min C ð PÞ

e ðm2 Þ is given by (5.4) and (5.16). We substitute C e ðm2 Þ into where C e k Þ to the form (5.17) and rearrange the problem ð P

Z

where

Z

e k ðTÞdx U

ð5:19Þ

X

ð5:7Þ

X

ð5:18Þ

X

e S k ðTÞ ¼



 ðdN : NÞdx þ e S k ðTÞ j ðN; TÞ 2 R ;

and



The tensor h in (5.6) is determined by the layout vY2 within Y. To find these compliances one can use the formula (3.4), which directly  gives the product: T a hab T b . One can prove that the formulation (5.4) is a full relaxation of (4.11).

e k ðTÞ ¼ minfhðhÞkTk2 þ kh j h 2 ½0; 1g: U

ð5:20Þ

e k ðTÞ is the convex envelope of the function U k ðTÞ One can note that U given by (4.16). We introduce the function

4872

S. Czarnecki et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4866–4881

Fig. 6. Transverse shear of the strip-wise plate.

~ k ðzÞ ¼ minfhðhÞz2 þ kh j h 2 ½0; 1g u

ð5:21Þ

and write (5.19) in the form

e S k ðTÞ ¼

Z

~ k ðkTkÞdx: u

ð5:22Þ

X

We rewrite (5.18) as below

min

Z

 ~ k ðkTkÞdx j ðN; TÞ 2 R : ½dN : N þ u

ð5:23Þ

By denoting

1 1 ~k ðkTkÞ; dN : N þ u 2 2

ð5:24Þ

we rearrange the problem (5.23) as follows:

min

Z X

 W k ðN; TÞdx j ðN; TÞ 2 R :

ð5:25Þ

If k is treated as fixed the problem above can be viewed as the equilibrium problem of an hyperelastic plate of a potential W k ðN; TÞ. We define the Lagrangian function

J ¼

Z X

fW k ðN; TÞ þ pv  Nab eab ðwÞ  T a ba ðw; vÞgdx;

ð5:26Þ

where v, w ¼ ðw1 ; w2 Þ play two roles: as test functions in (2.11) and Lagrangian multipliers of the statical admissibility. The stationary condition dJ ¼ 0 is satisfied for arbitrary dNab , dT a , dv, dwa . Due to the fields dNab , dT a being arbitrary, one finds

eab ðwÞ ¼

oW k ðN; TÞ oN

ab

;

ba ðw; vÞ ¼

oW k ðN; TÞ ; oT a

rffiffiffiffiffiffiffi k za ¼ Ha ; DH

ð5:27aÞ

a ¼ 1; 2; DH ¼ H2  H1

rffiffiffiffiffiffiffi k ; DH

b ¼ k

H1 : DH

oW k ðM; Q Þ oM

ab

;

ca ¼

oW k ðM; Q Þ ; oQ a

U 1 ¼ h1 z21 ;

U 2 ¼ h2 z22 þ k:

~ k ðkQ kÞ 1 ou : 2 oQ a

8 if kQ k 2 ð0; z1 Þ; h1 Q a > > < qffiffiffiffiffi a 1 k ca ¼ Q kQ k if kQ k 2 ðz1 ; z2 Þ; DH > > : a if kQ k P z2 : h2 Q

ð5:33Þ

The relation kckðkQ kÞ is shown in Fig. 8. The Eq. (5.33) can be inverted:

ð5:27bÞ

ð5:28Þ

~ k ðzÞ is the convex envelope of two parabolas (5.21), The function u see Fig. 7. It has the form:

8 2 if z 2 ð0; z1 Þ; > < h1 z ~ k ðzÞ ¼ az þ b u if z 2 ðz1 ; z2 Þ; > : h2 z2 þ k if z P z2 ;

ð5:32Þ

b k Þ of Section 4 We see now the relations between the problems ð P e k ðTÞ is the convexification of U k ðTÞ. In and (5.18). The function U the range: kTk 2 ð0; z1 Þ the core assumes the weakest properties, is made from the weaker material 1. If kTk > z2 the core assumes the strongest properties and is made from material 2. If kTk assumes the intermediate properties the core assumes the strip-wise composite microstructure, Fig. 6. Due to the potential W k ðM; Q Þ being positive and convex the problem (5.25) is well posed; its solution exists. Let us write the constitutive Eq. (5.28) explicitly in the subsequent regimes. We find

where qab ¼ eab ðuÞ, ca ¼ ba ðu; wÞ, see (2.7). The above Eq. (5.27b) can be written according to (5.24) as below

qab ¼ dabkl Mkl ; ca ¼

ð5:31Þ

One can prove that the straight line az þ b is tangent to both the parabolas: h1 z2 and h2 z2 þ k at points ðz1 ; U 1 Þ, ðz2 ; U 2 Þ where

while making use of the fields dv, dwa being arbitrary results in the equilibrium Eq. (2.11). The solution to the problem above will be denoted by N ¼ M, T ¼ Q , v ¼ w, w ¼ u. Let us re-write the stationarity conditions (5.27a) in the form

qab ¼

ð5:30Þ

and

a¼2

X

W k ðN; TÞ ¼

where

ð5:29Þ Fig. 7. Maxwell line for two parabolas.

4873

S. Czarnecki et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4866–4881

Fig. 8. The norm kck versus kQ k.

Fig. 10. Interpretation of the microstructure parameters.

aH ðu; w; w; vÞ ¼ f ðvÞ;

ð6:3Þ

valid for each kinematically admissible ðw; vÞ, where Z aH ðu; w; w; vÞ ¼ ½eab ðuÞDabkl ekl ðwÞ þ ba ðu; wÞHak ðm2 ; /Þbk ðw; vÞdx; X

ð6:4Þ

Fig. 9. The norm kQ k versus kck.

8 > h1 ca > > > < a Q ¼ fca > > > > :h c 2 a

qffiffiffiffiffi k ; DH qffiffiffiffiffi k if kck ¼ DH; qffiffiffiffiffi if kck > DkH:; if kck <

ð5:34Þ

e ¼ f ðwÞ is while f ðvÞ is defined by (2.12). The plate compliance C determined by the distribution of the design variables: m2 ðxÞ and /ðxÞ in domain X and not only by the distribution of m2 ðxÞ as in e problem ð PÞ. The considered problem of minimum compliance has the form

minff ðwÞjthe fieldsðw; u1 ; u2 Þ are solutions to ð6:3Þ and

where f satisfies H1 6 f 6 H2 . The relation kQ kðkckÞ is shown in Fig. 9.

At least two numerical approaches are possible. One can fix the multiplier k and elaborate a method for solving the non-linear ~ k ðÞ equilibrium problem (5.18), using the explicit expression of u given by (5.29). The value of m2 is fixed by the rule

¼ A2 ; 0 6 m2 6 1g ð6:5Þ

Lðu; w; w; v; m2 ; /; zþ ; z ; kÞ ¼ f ðwÞ þ f ðvÞ  aH ðu; w; w; vÞ Z  m2 dx  A2 þk X Z þ ½zþ ðm2  1Þ  z m2 dx;

ð6:1Þ

where zþ ðxÞ, z ðxÞ, k 2 R are Lagrangian multipliers. Here, v, w1 ,w2 play also the role of such multipliers, although they are also test functions in (6.3). The condition dL ¼ 0 with respect to dw, dua , dv, dwa , dm2 , d/, dzþ , dz , dk gives (a) The initial equations, namely Eq. (6.3) and the isoperimetric condition. (b) Necessary optimality conditions:

aH ðw; v; du; dwÞ ¼ f ðdwÞ; oHab ca c ¼ k þ z þ  z ; om2 b zþ P 0; zþ ðm2  1Þ ¼ 0;  z m2 ¼ 0; z P 0;

2

H11 ¼ H11 cos2 /  H12 sin 2/ þ H22 sin /; H

¼ ðH

11

H

22

Þ cos / sin / þ H

12

cos 2/;

ð6:6Þ

X

This method requires programming the equilibrium problem (5.18) or applying the user’s procedures, available e.g. in ABAQUS, for introduction of the own constitutive equations. The other possible approach is based on the interpretation of the underlying microstructure which realizes minimum in (5.8). Just this approach will be used in the sequel. The microstructure of Fig. 4 is determined by two parameters: the angle /ðxÞ of inclination of the strips and the area fraction m2 ¼ m2 ðxÞ, x 2 X (see Fig. 10). The transverse stiffnesses measured with respect to the basis ðe01 ; e02 Þ are given by the Eq. (3.9), where h2 ¼ m2 , h1 ¼ m1 . The components of stiffnesses with respect to the basis ðe1 ; e2 Þ are

12

m2 ðxÞdx

X

and the regularity assumptions are omitted. The Lagrangian functional for problem (6.5) reads

6. Numerical treatment of problem (5.17)

8 0 if kQ k 2 ð0; z1 Þ; > > <  qffiffiffiffiffi  1 D H m2 ¼ DH kQ k k  H1 if kQ k 2 ðz1 ; z2 Þ; > > : 1 if kQ k P z2 :

Z

ð6:2Þ

2

H22 ¼ H11 sin / þ H12 sin 2/ þ H22 cos2 /; while here H12 ¼ 0. These components are treated as functions Hab ðm2 ; /Þ where m2 ¼ m2 ðxÞ, / ¼ /ðxÞ. The bending stiffnesses are uniform in X. Substitution of (2.9) into (2.11) gives the variational equation

ca

oHab c ¼ 0: o/ b

ð6:7Þ ð6:8Þ ð6:9Þ ð6:10Þ ð6:11Þ

The fields dw, du are subject to the same kinematic conditions as w and u, hence (6.7) implies

w ¼ u;

v ¼ w;

ð6:12Þ

which is the known property of the minimum compliance problem. The condition (6.11) determines the optimal value of the angle

4874

S. Czarnecki et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4866–4881

n

~ / ~þ / 2 /;

po 2

ð6:13Þ

;

where

~¼ /

8 < 1 arctan 2 ðc

2c1 c2 2 ðc2 Þ2

1Þ

:p

if ðc1 Þ2 –ðc2 Þ2 ; if ðc1 Þ2 ¼ ðc2 Þ2 :

4

ð6:14Þ

We note that the formula (6.14) does not involve m2 and material constants. If 0 < m2 < 1 then zþ ¼ 0, z ¼ 0 and the condition (6.8) reduces to

ca

oHab c ¼k om2 b

ð6:15Þ

in which the derivatives oHab =om2 can be easy found analytically. In order to find optimal distributions of m2 and / we apply the COC method, cf. Bendsøe [12], applied already in the early papers by Cheng and Olhoff. The material of the faces will be assumed as isotropic, hence

e abkl ¼ C f



Ef 2

1  ðmf Þ

mf dab dkl þ

 1  mf ak bl ðd d þ dal dbk Þ ; 2

Step 4. Computing the stiffnesses Hab ¼ Hab ðm2 ; /Þ by (6.2), where Hab are given by (3.9) for ha ¼ ma . For k > 2 the values of Hab are computed for m2 , / found in Step 9. Step 5. The data of Steps 1 and 4 are stored and transmitted to the solver of statics. This solver is then executed. Step 6. The statics solver finds the components of strains ðqab Þ, ðca Þ and they are stored into a text file. We compute the total compliance. Step 7. Reading the data from the statics solver. By formula (6.14) ~ at each element. Of possible valwe find new values of / ues of /, see (6.13), we choose this value for which the energy density is the greatest. Step 8. For given values of / and m2 we compute values of the derivatives oHab =om2 . Step 9. We update m2 and k. Let for the index k

1 k

oHab ðm2 ; /ðkÞ Þ ðkÞ cb om2 ðkÞ

CðmðkÞ cðkÞ a 2 Þ ¼ and

ð6:16Þ

where Ef , mf are Young modulus and Poisson ratio of the material of the faces. The core is filled up by the honeycomb structure (with uniform dimensions l, t) made from two different isotropic materials of moduli ðE1 ; m1 Þ and ðE2 ; m2 Þ. These data determine the moduli la according to 2la ¼ Ea =ð1 þ ma Þ; the moduli lðcaÞ are given by (2.6) and the transverse stiffnesses – by (4.3). The equilibrium problem of the plate will be solved by the FEM with the triangular DSG3 element by Bletzinger et al. [33]. A critical analysis of this element has been performed recently in Kursa and Lewin´ski [35]. It turns out that this element is consistent. This means that for the regular meshes it produces the finite difference equations which converge to the partial differential equations of this plate theory, if the dimensions of the triangles tend to zero. Many popular triangular finite elements for moderately thick plates do not pass this self-evident consistency test, and this distinguishes this element from other 9 degrees of freedom triangular elements. This test, along with conventional tests based on direct comparisons with known analytical results for rectangular plates convince the present authors that this element is reliable and can be applied for the purpose of the minimum compliance design. In the DSG3 element the unknowns are: deflection w and two angles of rotations u1 , u2 at each node of the triangle. The bending deformations ðqab Þ and the transverse shear deformations ðca Þ are constant in the element. By analogy, the design variables m2 and / will be treated as constant within the element. Thus the approximation of the design variables will be element-wise constant, hence non-differentiable. Under these approximations concepts the COC algorithm was developed. It consists of the following steps:

sðkÞ ¼ m2 ðCðm2 ÞÞg : ðkÞ

ðkÞ

We find the new value of m2 by the scheme of Bendsøe [12]

ðkþ1Þ

m2

8 ðkÞ > > < maxfð1  nÞm2 ; 0g if ðiÞ; ðkÞ g ¼ mðkÞ if ðiiÞ; 2 ðCðm2 ÞÞ > > : ðkÞ minfð1 þ nÞm2 ; 1g if ðiiiÞ;

Fig. 11. Simply supported plate. Distribution of m2 .

Step 1. (a)

Assumption of the data: dimensions Lx , Ly of the rectangular plate, thicknesses c and d of the layers. Assumption of the FE triangulation of the domain: X ¼ ð0; Lx Þ  ð0; Ly Þ. The triangles are denoted by Xe , e ¼ 1; . . . ; N. Assumption of the transverse loading of intensity pðxÞ and its FE approximation. Assumption of the moduli Ef , mf of the faces and: ðE1 ; m1 Þ, ðE2 ; m2 Þ of the materials from which the cellular core is constructed. Fixing the kinematic boundary conditions. (b) Data for the COC method: initial values of m2 and / within all elements, value of k, error tolerance d, the move limit n and the damping factor g. Step 2. Computing the quantities common for the whole plate and all loops of the algorithm, such as the stiffnesses Dabkl . Step 3. Start the optimization loop, place the loop index k ¼ 1.

Fig. 12. Simply supported plate. Vector field ^c.

S. Czarnecki et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4866–4881

4875

Fig. 13. Hingely supported plate. Distribution of density m2 . Fig. 16. The same plate. Vector field ^c.

Fig. 14. Hingely supported plate. Vector field ^c.

Fig. 17. The plate clamped along two adjacent sides, subject to a uniform loading. Distribution of density m2 .

Fig. 15. The optimal clamped plate subject to a uniform loading. Distribution of density m2 .

Fig. 18. The same plate. Vector field ^c.

where ðkÞ

ðiÞ

sðkÞ 6 maxfð1  nÞm2 ; 0g;

ðiiÞ

maxfð1  nÞm2 ; 0g 6 sðkÞ 6 minfð1 þ nÞm2 ; 1g;

ðiiiÞ minfð1 þ

ðkÞ

ðkÞ nÞm2 ; 1g

ðkÞ

6 sðkÞ :

monotonic in k, which assures the convergence of the bisection method. Step 10. We store new values of m2 at each elements, change the index k into k þ 1 and go to Step 4.

ðkþ1Þ m2

Upon computing at each element we compute the integral R m dx according to the assumed element-wise constant distribu2 X tion of m2 . By the bisection method we find new multiplier k which R assures the isoperimetric condition X m2 dx ¼ A2 to be fulfilled with given accuracy d. We take advantage of the integral of m2 being

7. Exemplary optimal designs of the core layer We consider a rectangular plate of dimensions Lx  Lx , Lx ¼ 500 cm. The thicknesses of the layers are: c ¼ 16 cm and

4876

S. Czarnecki et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4866–4881

Fig. 19. The plate clamped along two adjacent sides, subject to a point load at the free corner. Distribution of density m2 . Fig. 20. The same plate. Vector field ^c.

d ¼ 0:3 cm; h ¼ 32:6 cm. The faces are made from steel of moduli Ef ¼ 210 GPa and mf ¼ 0:21. The core has the honeycomb structure of regular hexagonal cells of the vertical walls of length l ¼ 1 cm, height 2c ¼ 32 cm and thickness t ¼ 0:02 cm. We assume that the honeycomb structure is made from two different isotropic materials: the first material is aluminum ðE1 ¼ 69:64 GPa;m1 ¼ 0:31Þ and the second is glass/epoxy, also treated as isotropic ðE2 ¼ 38:96 GPa;m2 ¼ 0:28Þ. The data above determine the moduli ð2Þ l1 and l2 and then the moduli lð1Þ c , lc by (2.6) and the transverse stiffness tensors H1 and H2 by (2.5). At the data assumed the con-

trast of stiffness equals: H2 =H1 ¼ 1:746. It is assumed that the material 2 occupies 51% of the domain or A2 =jXj ¼ 0:51. Example A. We consider a simply supported plate:

w ¼ 0;

M n ¼ 0;

Ms ¼ 0

ð7:1Þ 2

along the edges. The loading is uniform: p ¼ 0:001 kN=cm . After 100 iterations one finds the distribution of m2 within the core layer, as in Fig. 11.

Fig. 21. The plate lying on four square supports, subject to a uniform load. Distribution of density m2 is shown in (b) and (c). Types of materials are shown in (d): stronger material (black) and first-rank laminate (gray).

S. Czarnecki et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4866–4881

4877

Example B. Consider a hingely supported plate, with other data taken as in Example A. The one quarter of the plate was divided into 30 by 30 elements, hatched region of the plate, see Fig. 13a. The boundary conditions are:

w ¼ 0;

M n ¼ 0;

us ¼ 0:

ð7:2Þ

This support is stiffer than the simple support, both the supports differing at most in the vicinity of the corners. Note, that the Kirchhoff plate theory does not distinguish these two kinds of boundary conditions. The hinge support at the corners is so stiff that the softer material is placed there, see Fig. 13. The strip directions are now different mainly in the vicinity of the corners, cf. Figs. 13 and 14.

Fig. 22. The same plate. Vector field ^c.

The computations were performed at various triangulations of the one quarter of the plate, obtaining similar results. The results of Figs. 11 and 12 were obtained with using 30 by 30 triangular DSG3 elements. The transverse shear deformation vector ^c ¼ ðc2 ; c1 Þ forms the vector field which indicates the directions orthogonal to the strips of the underlying microstructure, see Figs. 12 and 6. One notes that the strips go along the diagonals of the plate.

Example C. Consider the plate clamped along all edges, see Fig. 15, subject to a uniform loading. One quarter of the plate was divided into 1800 elements. We note that the zones around the corners are occupied by the weaker material. The vector field ^c has similar distribution as in Example B, see Fig. 16. Example D. Consider the plate clamped along two adjacent edges and free on the two remaining edges, see Fig. 17. Along the clamped edges we have:

w ¼ 0;

un ¼ 0; us ¼ 0:

ð7:3Þ

The whole plate was divided into 50 by 50 triangular elements. The loading is the same as before. We note that the stronger material is

Fig. 23. The plate lying on four square supports at the corners, subject to a uniform load: (a) the static scheme and the one quarter considered; (b), (c) optimal layouts of m2 in the whole plate and in one quarter and (d) types of materials – stronger material (black) and first-rank laminate – gray.

4878

S. Czarnecki et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4866–4881

one quarter was divided into 60 by 60 elements. The strips around the edges of the supports occur to be orthogonal to their square contours. The stronger material is located around the supports, see Figs. 21 and 22.

Fig. 24. The same plate. Vector field ^c.

not necessary in the zone around the clamped corner and along the diagonal. The strips are orthogonal to the clamped edges and parallel to the free edges, cf Fig. 18. Example E. We consider the plate supported as in Example D, subject to a point load at the free corner. The triangulation was the same as in the previous example (see Fig. 18). Optimization stiffens the plate along the free edges, see Fig. 19. The strips go along these edges, see Fig. 20. Example F. The plate lies at four square-shaped supports of dimensions 25  25 cm, their centers being at distance 87.5 cm to the nearest edge. The nodes within the square supports are fixed. The loading is uniform of intensity p ¼ 0:001 kN=cm2 . The

Example G. The plate lies at four square-shaped supports of dimensions 25  25 cm, their centers being at near the corners. The loading is uniform of intensity p ¼ 0:001 kN=cm2 . The triangulation was the same as in the previous example. The strips around the edges of the supports are orthogonal to their square contours. The stronger material is located around the supports, see Figs. 23 and 24. In all the examples the final layout was drastically different than the initial one, but this essential change in the layout was not accompanied by the essential change of the compliance, measured with respect to the compliance of the initial homogeneous design. For instance, the relative decay of the compliance for the design in Fig. 15 equals 4%, while the layout in Fig. 21 is associated with the decay of 9.3%. Example H. Consider the plate simply supported along all the edges with the boundary conditions (7.1). The plate is subject to four point loads at the centers of the plate quarters, as shown in Fig. 25a. The material and geometrical data are the same as in the previous examples. One quarter of the plate was divided into 5000 elements. We note that the zones around the point loads are occupied by the stronger material. The vector field ^c circulates around the point loads, see Fig. 26. The problem of optimal design of Reissner–Mindlin plates of piece-wise varying thickness was discussed in Ref. [22]. As a result of relaxation of the variational formulation of the optimal design

Fig. 25. Simply supported plate subject to the four point loads: (a) the static scheme and the one quarter considered; (b), (c) optimal layouts of m2 in the whole plate and in one quarter and (d) types of materials – stronger material (black) and first-rank laminate – gray.

S. Czarnecki et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4866–4881

Fig. 26. The same plate. Vector field ^c.

problem of the compliance minimization the composite domains occur of the underlying laminated microstructure of first-, secondand third-rank. The moment formulation was applied and the theoretical consideration was illustrated by an example of a simply supported plate with the same loading as in Fig. 25a. The design variables affected both the bending and transverse shear stiffnesses. Therefore, the optimal layout of the thicknesses given in Figs. 6a, 7a in [22] assumes the form of the diagonal cross-like structure which links the opposite corners of the plate. This solution results mainly from adaptation of the plate characteristics to the shape of the moment diagrams in order to stiffen the domains where the bending deformation prevails. One can suspect that the deflection due to transverse shear are of minor importance. In Ref. [22] two cases of optimal control were analyzed in which the laminates

4879

up to third-rank and then – up to second-rank were admitted, see Figs. 6b, 7b in [22]. It turned out that the domain where the first-rank laminate occurred is almost insensitive to the choice mentioned. The first-rank laminate surrounds the domains where the maximal value of the thickness is achieved. Due to the loading and the boundary conditions being the same in both the examples: of the sandwich plate of constant thickness and the composite core considered in the present paper (Figs. 25 and 26) and the plate of piece-wise varying thickness considered in Ref. [22] it is thought appropriate to compare the final optimal solutions. The core of the sandwich plate considered here adapts to the distribution of the transverse shear forces while the distribution of the moments does not affect directly the optimal layout. The reason is that the transverse forces in the plate of Fig. 25a are negligible in the central domain of the plate and in the vicinity of the corners. The optimal distribution of m2 assumes the maximal values in the areas around the point loads and transmits them to the edges along certain arch-like regions, see Fig. 25d, being the domains in which the function m2 achieves its maximal value 1. It is worth comparing the distribution of m2 of Fig. 25b with the vector field graph ^c from Fig. 26 to note that in the regions where m2 is almost zero the values of the transverse shear deformations are minimal. Therefore, the layout of m2 of Fig. 6a in [22] is inevitably different than the layout of m2 of Fig. 25b of the present paper. Example I. Consider the plate hingely supported along all the edges with the boundary conditions (7.2). The plate is subject to four point loads at the centers of the plate quarters, as shown in Fig. 27a. The material and geometrical data are the same as in the previous examples. The vector field ^c circulates around the point loads, see Fig. 28.

Fig. 27. Hingely supported plate subject to the four point loads: (a) the static scheme and the one quarter considered; (b), (c) optimal layouts of m2 in the whole plate and in one quarter and (d) types of materials – stronger material (black) and first-rank laminate – gray.

4880

S. Czarnecki et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4866–4881

[2] [3] [4] [5]

[6] [7] [8]

Fig. 28. The same plate. Vector field ^c.

The boundary conditions of a hinge support (7.2) are stiffer than the boundary conditions of a simple support (7.1) and that is why the corner domains of the plate of Fig. 27a do not need stiffening, which is reflected by the m2 values achieving the minimal values. The stronger material (m2 ¼ 1) is located so as to transmit the point loads to the edges along the arch-like bands, see Fig. 27d. The vector field ^c in Fig. 28 is different than in the problem of the simply supported plate, cf. Fig. 26. The main discrepancies occur in the corner domains.

[9]

[10]

[11] [12] [13] [14] [15]

8. Final remarks Due to a mathematical analogy between the second local homogenization problem of the plates of moderate thickness and the well known local homogenization problem of the model elliptic scalar problem the relaxation of problem ðPÞ could be performed by the known methods of topology optimization. Note, however, that the general problem of relaxation of the optimal transversely homogeneous layout of two materials in the moderately thick Reissner–Hencky plates has not been solved till now. The reason lies in the difficulties arising from the bending-transverse shearing coupling; the design variable v2 appears in both the tensors D and H, which makes it impossible to decouple the problem. In the problem considered, in which the two-component core layer was designed, the first-rank composites occurred to be sufficient to saturate all the bounds and make the minimum of the compliance achievable. This problem is relatively easier than the problem of designing the whole plate in which two stiffness tensors are subject to optimization. Thank of it the final results are less sophisticated than those known from plane elasticity in which at least second-rank composite should be a taken into account. Consequently, the results presented are also easier to manufacture, since the strip-wise domains can be approximated by introduction of appropriate ribs or stiffeners. Acknowledgements The work was partially supported within the grant 4T07A 038 30 of the Polish Ministry of Science and Higher Education: Theory and numerical implementation of relaxed formulations of optimization problems with coupled fields. Designing of the layout of materials in the composite structures. References [1] L.V. Gibiansky, A.V. Cherkaev, Microstructures of elastic composites of extremal stiffness and exact estimates of the energy stored in them, in: A.V.

[16]

[17]

[18]

[19]

[20] [21]

[22]

[23] [24] [25] [26] [27] [28] [29] [30]

[31]

Cherkaev, R.V. Kohn (Eds.), Fiziko-Tekhnichesk. Inst. Im. A.F. Ioffe. AN SSSR, Preprint No. 1115 (1987). Leningrad (in Russian), p. 52. English Translation, Topics in the Mathematical Modelling of Composite Materials, Birkhäuser, Boston, 1997. A. Cherkaev, Variational Methods for Structural Optimization, Springer, New York, 2000. G. Allaire, Shape Optimization by the Homogenization Method, Springer, New York, 2002. R. Lipton, On a saddle-point theorem with application to structural optimization, J. Optim. Theory Appl. 81 (1994) 549–568. J.B. Jacobsen, N. Olhoff, E. Rønholt, Generalized shape optimization of threedimensional structures using materials with optimum microstructures, Mech. Mater. 28 (1998) 207–225. T. Borrvall, J. Petersson, Large-scale topology optimization in 3D using parallel computing, Comput. Methods Appl. Mech. Engrg. 190 (2001) 6201–6229. A. Dı´az, R. Lipton, Optimal material layout for 3D elastic structures, Struct. Optim. 13 (1997) 60–64. N. Olhoff, E. Rønholt, J. Scheel, Topology optimization of threedimensional structures using optimum microstructures, Struct. Optim. 16 (1998) 1–18. S. Czarnecki, T. Lewin´ski, Shaping the stiffest three-dimensional structures from two given isotropic materials, Comput. Assist. Mech. Engrg. Sci. 13 (2006) 53–83. G. Allaire, R.V. Kohn, Optimal design for minimum weight and compliance in plane stress using extremal microstructures, Eur. J. Mech., A/Solids 12 (1993) 839–878. G. Allaire, E. Bonnetier, G. Francfort, F. Jouve, Shape optimization by the homogenization method, Numer. Math. 76 (1997) 27–68. M.P. Bendsøe, Optimization of Structural Topology, Shape, and Material, Springer, New York, 1995. M.P. Bendsøe, O. Sigmund, Material interpolation schemes in topology optimization, Arch. Appl. Mech. 69 (1999) 635–654. M.P. Bendsøe, O. Sigmund, Topology Optimization, Theory, Methods and Applications, Springer, Berlin, 2003. L.V. Gibiansky, A.V. Cherkaev, Designing composite plates of extremal rigidity, in: A.V. Cherkaev, R.V. Kohn (Eds.), Fiziko-Tekhnichesk. Inst. Im. A.F. Ioffe. AN SSSR, Preprint No. 914 (1984). Leningrad (in Russian). English Translation, Topics in the Mathematical Modelling of Composite Materials, Birkhäuser, Boston, 1997. T. Lewin´ski, J.J. Telega, Plates, Laminates and Shells, Asymptotic Analysis and Homogenization, World Scientific, Singapore, New Jersey, London, Hong Kong, 2000. L.A. Krog, N. Olhoff, Topology and reinforcement layout optimization of disk, plate, and shell structures, in: G.I.N. Rozvany (Ed.), Topology Optimization in Structural Mechanics, Springer, Wien-New York, 1997, pp. 237–322. S. Czarnecki, T. Lewin´ski, Optimal layouts of a two-phase isotropic material in thin elastic plates, in: Z. Waszczyszyn, J. Pamin, (Eds.), Proceedings of the Second European Conference on Computational Mechanics, ECCM-2001, Cracow, 26-29 VI 2001, CD ROM, 2001. K. Kolanek, T. Lewin´ski, Circular and annular two-phase plates of minimal compliance, Comput. Assist. Mech. Engrg. Sci. 10 (2003) 177– 199. _ G. Dzierzanowski, Optimization of the layout of materials in elastic plates and shells, Ph.D. Thesis. Warsaw, 2004 (in Polish). T. Lewin´ski, Homogenization and optimal design in structural mechanics, in: P.P. Castañeda, J.J. Telega, B. Gambin (Eds.), Nonlinear Homogenization and its Application to Composites, Polycrystals and Smart Materials, NATO Science Series II, Mathematics, Physics and Chemistry, vol. 170, Kluwer, Dordrecht, 2004, pp. 139–168. A.R. Dı´az, R. Lipton, C.A. Soto, A new formulation of the problem of optimum reinforcement of Reissner–Mindlin plates, Comput. Methods Appl. Mech. Engrg. 123 (1995) 121–139. G.A. Francfort, F. Murat, L. Tartar, Fourth-order moments of nonnegative measures on and applications, Arch. Rat. Mech. Anal. 131 (1995) 305–333. E. Reissner, On bending of elastic plates, Quart. Appl. Math. 5 (1947) 55–68. T. Lewin´ski, On displacement-based theories of sandwich plates with soft core, J. Engrg. Math. 25 (1991) 223–241. N. Buannic, P. Cartraud, T. Quesnel, Homogenization of corrugated core sandwich panels, Compos. Struct. 59 (2003) 299–312. L.J. Gibson, M.F. Ashby, Cellular Solids, second ed., Structure and Properties, Cambridge University Press, 1999. G. Shi, P. Tong, Equivalent transverse shear stiffness of honeycomb cores, Int. J. Solids Struct. 32 (1995) 1383–1393. W. Becker, Closed-form analysis of the thickness effect of regular honeycomb core material, Compos. Struct. 48 (2000) 67–70. J. Hohe, A direct homogenization approach for determination of the stiffness matrix for microheterogeneous plates with application to sandwich panels, Compos.: Part B 34 (2003) 615–626. A. Bourgeat, R. Tapiéro, Homogénéisation d’une plaque mince, thermoélastique, perforée transversalement, de structure non uniformement périodique, dans le modèle de la théorie naturelle, C.R. Acad. Sci. Paris. Sér. I 297 (1983) 213–216.

S. Czarnecki et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 4866–4881 [32] T. Lewin´ski, Homogenizing stiffnesses of plates with periodic structure, Int. J. Solids Struct. 29 (1992) 309–326. [33] K.U. Bletzinger, M. Bischoff, E. Ramm, A unified approach for shear-lockingfree triangular and rectangular shell finite elements, Comput. Struct. 75 (2000) 321–334. [34] J. Rakowski, A new methodology of evaluation of C 0 bending finite elements, Comput. Methods Appl. Mech. Engrg. 91 (1991) 1327–1338.

4881

[35] M. Kursa, T. Lewin´ski, Evaluation of convergence of DST and DSG3 elements in equilibrium problems of elastic plates of moderate thickness, in: K. Dems, Z. Wie ß ckowski, (Eds.), Proceedings of the 17th International Conference on Computer Methods in Mechanics (CMM 2007) Spała, Poland, June 19–22, 2007. [36] J. Necˇas, M. Hlavacˇek, Mathematical Theory of Elastic and Elasto-Plastic Bodies, An Introduction, Elsevier, Rotterdam, 1981.