Topology optimization of incompressible materials based on the mixed SBFEM

Computers and Structures 165 (2016) 24–33

Contents lists available at ScienceDirect

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

Topology optimization of incompressible materials based on the mixed SBFEM Chao Li, Liyong Tong ⇑ School of Aerospace, Mechanical and Mechatronic Engineering, The University of Sydney, NSW 2006, Australia

a r t i c l e

i n f o

Article history: Received 27 July 2015 Accepted 1 December 2015

Keywords: Topology optimization Incompressible materials Mixed scaled boundary finite element method

a b s t r a c t This paper presents the topology optimization based on a semi-analytical mixed scaled boundary finite element method (SBFEM) for problems involving incompressible materials using a moving iso-surface threshold method (MIST). The mixed SBFEM does not have volumetric issues in the presence of incompressibility and it can pass the inf–sup stability conditions without requiring any mesh-dependent and stabilization parameters. Thus the topology optimization of incompressible materials is also possible based on the mixed SBFEM. The solutions of stress and strain in the mixed SBFEM are derived semianalytically. A response function in MIST is chosen as the mutual strain energy density which is formulated by the stress and the strain obtained in the mixed SBFEM. Numerical examples including the optimal design problems of compressible and incompressible materials as well as the rubber mounts containing incompressible solid and fluid are presented. The present results are also verified with the available reference results. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Incompressible and nearly incompressible behavior are often encountered in the practical engineering problems. For example, most elastomers and rubber like materials have negligible compressible effects. It is well known that although the standard displacement-based finite element method (FEM) is widely used for solving linear and non-linear elastic problems, it exhibits volumetric locking issues when Poisson’s ratio approaches to 0.5. Thus, it is also difficult to conduct the standard finite element based topology optimization of the incompressible materials. Since Herrmann [1] proposed a mixed FEM formulations to overcome volumetric locking issues, many techniques such the selective reduced integration method [2], the B-bar method [3] and the assumed strain methods [4–6] have been developed. However, the standard mixed formulations [7,8] cannot produce stable results of pressure field near the incompressible limit, which is referred to the Ladyzhenskaya–Babuška–Brezzi (LBB) inf–sup stability conditions [9,10]. Later, the stabilized methods with the use of stabilization parameters was developed to avoid the drawbacks in the mixed FEM [11,12]. A mixed formulation based on the semi-analytical SBFEM [13] was developed in Li and Tong [14] to solve the fracture problems of incompressible materials. The mixed SBFEM will be used to investigate the optimal design ⇑ Corresponding author. Tel.: +61 2 93516949; fax: +61 2 93514841. E-mail address: [email protected] (L. Tong). http://dx.doi.org/10.1016/j.compstruc.2015.12.003 0045-7949/Ó 2015 Elsevier Ltd. All rights reserved.

of the structures involving the incompressible materials in this paper. Topology optimization of incompressible media using density approach (SIMP) has been substantially studied based on the mixed FEM. Sigmund and Clausen [15] applied a mixed u=p finite element formulation which uses an incompressible medium to transfer pressure loads to solve the pressure load problem using the density approach to topology optimization. Bruggi and Venini [16] proposed an alternative mixed formulation consisting in a truly-mixed variational formulation coupled to a mixed finite element formulation for the topology optimization of incompressible body. This formulation using the triangular elements of Johnson and Mercier [17] with the stresses as main variables can pass the inf–sup stability condition. Bruggi and Cinquini [18] also extended this formulation for the topology optimization involving fluid incompressibility considering pressure loads. Jang and Kim [19] employed the displacement-based nonconforming FEM for the topology optimization of structures involving incompressible materials by minimizing the mean compliance and Panganiban et al. [20] applied the displacement-based nonconforming FEM for the design of pressure-actuated compliant mechanisms. Jang et al. [21] used the displacement-based P1-nonconforming FEM for the optimization of pressure-loaded structure problems. Alternatively, the moving iso-surface threshold topology optimization method (MIST) developed by Tong and Lin [28] was used in conjunction with a mixed u=p finite element formulation to maximize structural stiffness for the optimal design of cellular

25

C. Li, L. Tong / Computers and Structures 165 (2016) 24–33

4

4

2

2

(a)

(b) F=2 Fig. 2. Design domain of a bridge-like structure subjected to one point load.

(c) Fig. 1. The scaled boundary finite element method: (a) typical mesh; (b) a subdomain with the scaling center O, radial coordinate n and local coordinate g; and (c) linear element on boundary.

pressurized morphing structures [22,23]. In MIST, the optimal topologies can be generated without requiring sensitivity analysis [24]. This method has also been applied to investigate shape morphing for piezoelectric structures [25]. In this paper, we extend the mixed SBFEM to solve the topology optimization problem of incompressible materials with the use of MIST method. The mixed SBFEM does not lock near the incompressible limit and it can pass the inf–sup stability conditions without requiring any mesh-dependent and stabilization parameters. Unlike in the FEM, the mixed SBFEM retains the features of semianalytical SBFEM which only discretizes the boundaries of the problem domain so that the modeled spatial dimensions are reduced by one. Standard finite element shape functions are used for both pressure and displacement fields on the boundary and they are interpolated identically. Moreover, the solutions of displacement, stress and strain are expressed analytically in the radial direction, which leads to more accurate results. Its accuracy has been demonstrated in Li and Tong [14]. In the MIST, a physical response function is calculated to update the design variables. In this study, the physical response function is chosen as the mutual strain energy density formulated by the semi-analytical solutions of the stress and the strain in the mixed SBFEM. This paper is organized as follows. The basic concept of the mixed SBFEM and its solution are first presented in Section 2. The MIST method is introduced in Section 3. Numerical studies including the widely studied optimal design problems for compressible and incompressible materials are presented in Section 5. Finally, the conclusions and discussions are addressed in Section 5. 2. The mixed SBFEM formulations In this section, the mixed SBFEM formulations are briefly stated. The detailed derivation is referred to Li and Tong [14].

Static equilibrium condition of the problem domain is enforced in the mixed form based on the virtual work principle, in which the internal virtual work equals to the external virtual work [26]

Z

Z

fdegT ½Dd fegdV þ V

fug ¼ ½ux ; uy T is the displacement vector. ft g is the surface traction. G is shear modulus which is given in terms of the Young’s modulus E and the Poisson’s ratio m and is independent on dimension with E G ¼ 2ð1þ mÞ. E is the Young’s modulus and m is Poisson’s ratio. The diagonal matrix ½I0  is

2 3 2 0 0 16 7 ½I 0  ¼ 4 0 2 0 5 2 0 0 1

ð3Þ

p is the hydrostatic pressure

p¼

1 fmgT frg 2

ð4Þ

and it is related to the volume strain ev by the bulk modulus K

p ¼ K ev ¼ KfmgT feg

ð5Þ

where fmg ¼ ½1; 1; 0T . The bulk modulus K depends on dimension and strain assumptions. For 3D, 2D plane strain and plane stress, E E the bulk modulus is given as K ¼ 3ð12 mÞ, K ¼ 2ð1þmÞð12mÞ and E K ¼ 2ð1 mÞ. The constitutive law in the mixed form for 2D linear

incompressibility is written as

frg ¼ ½Dd feg þ fmgp

ð6Þ

 T in which frg ¼ rxx ; ryy ; rxy denotes the stress vector. In linear elasticity, feg is given as

feg ¼ ½Lfug

ð7Þ

where ½L is the linear differential operator [13]. Without considering the body force, the equilibrium equation in the domain is written as


g ¼ 0 ½LT fr

ð8Þ

2.2. The mixed SBFEM formulations

2.1. Problem statement

Z

which is valid in the volume V bounded by the surface vector S. In h iT above equation, feg ¼ exx ; eyy ; cxy is the strain vector.

fdegT fmgfpgdV  V

fdugT ftgdg ¼ 0

ð1Þ

S

with

  1 ½Dd  ¼ 2G ½I0   fmgfmgT 2

ð2Þ

In the SBFEM, a problem domain is subdivided into multisubdomains as shown in Fig. 1(a). In Fig. 1(b), a scaling center O is chosen in the subdomain from which entire boundary is visible. A radial coordinate n is defined from the scaling center ðn ¼ 0Þ to the boundary ðn ¼ 1Þ. The nodal displacement function fuðnÞg and the pressure function fpðnÞg are introduced in the radial coordinate along the radial lines passing through the scaling center O and a node on the boundary. In this study, linear line-elements are used with the local coordinate g varying from 1 to +1 as shown in Fig. 1(c). The displacement fuðn; gÞg and pressure

26

C. Li, L. Tong / Computers and Structures 165 (2016) 24–33

where ½I is a 2  2 identity matrix. The strain field is related to the displacement field by [27]

h i i 1h feðn; gÞg ¼ B1 ðgÞ fuðnÞg;n þ B2 ðgÞ fuðnÞg n

ð12Þ

The stress at any point frðn; gÞg can be obtained as

 i i  h 1h B1 ðgÞ fuðnÞg;n þ B2 ðgÞ fuðnÞg frðn; gÞg ¼ ½Dd  þ K fmgfmgT n ð13Þ The pressure denoted as** fpðn; gÞg is analogously obtained as

fpðn; gÞg ¼ K fmgT

(a)

h

 i i 1h B1 ðgÞ fuðnÞg;n þ B2 ðgÞ fuðnÞg n

ð14Þ

Based on the virtual work principle, the governing equation of the mixed SBFEM is derived as T

½E0 n2 fuðnÞg;nn þ ð½E0   ½E1  þ ½E1  ÞnfuðnÞg;n  ½E2 fuðnÞg ¼ 0

ð15Þ

40*40 SBFE subdomains T

fF ðnÞg ¼ ½E0 nfuðnÞg;n þ ½E1  fuðnÞg

ð16Þ

with

Z

þ1

½E0  ¼ ½E1  ¼ ½E2  ¼

1 Z þ1 1 Z þ1 1

(b)

h iT  i h B1 ðgÞ ½Dd  þ fmgK fmgT B1 ðgÞ jJjdg

ð17aÞ

h iT  i h B2 ðgÞ ½Dd  þ fmgK fmgT B1 ðgÞ jJjdg

ð17bÞ

h iT  i h B2 ðgÞ ½Dd  þ fmgK fmgT B2 ðgÞ jJjdg

ð17cÞ

Apply block-diagonal Schur decomposition, resulting in the general solutions for the displacements and the internal nodal forces to the Eq. (15)

fuðnÞg ¼

N h X

i

WiðuÞ n½Si  fci g

ð18aÞ

i¼1

80*80 SBFE subdomains

fF ðnÞg ¼

N h X

i

WðqÞ n½Si  fci g i

ð18bÞ

i¼1

where ½Si  ði ¼ 1; 2; . . . ; NÞ are the diagonal blocks of the real Schur h i ðuÞ matrix and Wi are the corresponding displacement modes. fci g ði ¼ 1; 2; . . . ; NÞ are the integration constants, which are determined by solving Eq. (18a) on the boundary. Eliminating the integration constant fci g from Eqs. (18a) and (18b) on the boundary results in

½Kfug ¼ fFg

(c)

ð19Þ

where the static stiffness matrix ½K is equal to

h i h i1 ðuÞ ½K ¼ ½W1ðqÞ     ½WNðqÞ   ½WðuÞ    ½WN  1  Fig. 3. Objective function histories for compressible design using different meshes: (a) 20  20; (b) 40  40; and (c) 80  80 SBFEM subdomains.

fpðn; gÞg at any point ðn; gÞ inside a domain are obtained by interpolating fuðnÞg and fpðnÞg with the linear shape functions

fuðn; gÞg ¼ ½Nu ðgÞfuðnÞg ¼ ½N1 ðgÞ½I; N2 ðgÞ½IfuðnÞg fpðn; gÞg ¼ ½Np ðgÞfpðnÞg ¼ ½N 1 ðgÞ; N2 ðgÞfpðnÞg

The strains and stresses are evaluated at a specified local coordinate

g within an element by substituting Eq. (18a) into Eqs. (12) and (13) feðn; gÞg ¼

ð9Þ ð10Þ

1 ð1  gÞ 2 1 N2 ðgÞ ¼ ð1 þ gÞ 2

N1 X ½Wei ðgÞn½Si ½I fci g

ð21Þ

i¼1

frðn; gÞg ¼

N1  X

 ½Dd  þ K fmgfmgT ½Wei ðgÞn½Si ½I fci g

ð22Þ

i¼1

with

N1 ðgÞ ¼

ð20Þ

ð11aÞ ð11bÞ

with

  h i ðuÞ ðuÞ ½Wei ðgÞ ¼ ½B1 ðgÞ Wi ½Si  þ ½B2 ðgÞ½Wi 

ð23Þ

Substitute Eq. (18a) into Eqs. (12) and (14) leading to the pressure

27

C. Li, L. Tong / Computers and Structures 165 (2016) 24–33

Phi function

Phi function

iso-surface

iso-surface

(a)

(b)

Phi function

Phi function

iso-surface

iso-surface

(c)

(d)

Fig. 4. / functions and iso-surfaces for compressible design of one point loaded bridge-like structure at iterations: (a) 5; (b) 10; (c) 40; and (d) 60.

fpðn; gÞg ¼

N1 X K fmgT ½Wei ðgÞn½Si ½I fci g

ð24Þ

i¼1

K ðxe Þ ¼ K void þ ðxe Þb ðK solid  K void Þ

ð25aÞ

Gðxe Þ ¼ Gvoid þ ðxe Þb ðGsolid  Gvoid Þ

ð25bÞ

with 3. Formulation of topology optimization Topology optimization aims to find the distribution of solid materials by minimizing or maximizing an objective function. In the present study, MIST approach is used for the topology optimization of incompressible media [28,24]. 3.1. Material interpolation For the conventional compressible material, the standard way is using a SIMP interpolation on the Young’s modulus E. As proposed in Sigmund and Clausen [15], the interpolation scheme for incompressible material is based on the bulk modulus K and shear modulus G and they are interpolated as the functions of the element density design variables, i.e., the element weighting factor xe . For a two-phase incompressible solid/void, the following interpolation in each element is used

0 6 xe 6 1

ð26Þ

where b is the penalization factor. K solid /Gsolid and K void =Gvoid are the bulk modulus/the shear modulus of the solid material and void, respectively. In the present study, the value of K solid is chosen to be 100 times larger than Gsolid for incompressible material. In order to avoid the singular issue, K void and Gvoid are set to be 0.001. In this interpolation scheme, the material is represented as void when xe ¼ 0 and as solid when xe ¼ 1. When considering a two phase solid/fluid interpolation, the material is represented as a fluid when xe ¼ 0 and a solid when xe ¼ 1. The interpolation scheme is represented as

K ðxe Þ ¼ K solid ¼ K fluid  Gðxe Þ ¼ Gfluid þ ðxe Þb Gsolid  Gfluid

ð27aÞ ð27bÞ

28

C. Li, L. Tong / Computers and Structures 165 (2016) 24–33

Phi function

Phi function

iso-surface

iso-surface

(a)

(b)

Phi function

Phi function

iso-surface

iso-surface

(c)

(d)

Fig. 5. / functions and iso-surfaces for incompressible design of one point loaded bridge-like structure at iterations: (a) 5; (b) 10; (c) 40; and (d) 60.

4

2

4

F=2

2

F=2 Fig. 7. Design domain of a bridge-like structure subjected to two point loads.

3.2. MIST method

Fig. 6. Objective function histories for the single-point loaded bridge-like incompressible structure.

MIST method developed by Tong and Lin [28] is as an effective topology optimization approach. Material distribution of design domain is obtained by minimizing or maximizing an objective function. The design domain consists of a regular SBFEM mesh as

29

C. Li, L. Tong / Computers and Structures 165 (2016) 24–33

(a)

Fig. 10. Dynamic filter size.

(b) (a) Fig. 8. Objective function histories for the two-point loaded bridge-like structure: (a) compressible material and (b) incompressible material.

F=2

1

(b) F=2 1 Fig. 9. Design domain of a cantilever beam subjected to symmetric loads.

shown in Fig. 1(a). The material property of each subdomain is evaluated by the element weighting factor xe which is between 0 and 1. The elements with a weighting factor of 1 represent solid region whereas the elements with a weighting factor of 0 represent void or fluid region.

Fig. 11. Objective function histories for the cantilever beam: (a) compressible material and (b) incompressible material.

3.2.1. Objective function When designing structures, the stiffness or the compliance must be taken into account. The formulation of the compliance is expressed as:

C ¼ fF gT fug

ð28Þ

30

C. Li, L. Tong / Computers and Structures 165 (2016) 24–33

(a) fixed solid region 1

1 Fig. 12. Design domain of a fluid-rubber mount structure subjected to a uniform distributed pressure.

A typical compliant mechanism problem is to maximize the displacement of the output DOF ui with a material volume fraction constraint, which can be achieved by maximizing the total mutual strain energy. Therefore, the energy objective function can be expressed as

Z

max : uout ¼ V

T

frðn; gÞgout feðn; gÞgin dV

(b)

ð29aÞ

Subject to : ½K fug ¼ fF g R x dV RV e 6 f solid dV V

ð29bÞ

0 6 xe 6 1

ð29dÞ

ð29cÞ

where the subscripts ‘in’ and ‘out’ refer to input and output load cases, respectively. f solid is the volume fractions of the solid material. 3.2.2. Selection of the response function / In MIST, a physical response function / is calculated at each node across the design domain to find the optimal topology of the structure and to update the element weighting factors during the iterative process. The response function for this optimization problem is the mutual strain energy density expressed as:

/ ¼ frgTout fegin

ð30Þ

Substituting Eqs. (25) or (27), (21) and (22) into Eq. (30) leads to

!2 N1   X ½Wei ðgÞfci g / ¼ ½Dd ðxe Þ þ K ðxe ÞfmgfmgT i¼1

ð31Þ

Fig. 14. Objective function histories for the incompressible fluid-rubber mount structure: (a) volume fraction of solid material 45% and (b) volume fraction of solid material 60%.

3.3. MIST algorithm The topology optimization problems are solved using an inhouse code written in MATLAB R2013b on an Intel(R) Core(TM) i7 CPU @ 4.0 GHz desktop machine. The MIST algorithm based on the mixed SBFEM can be described as follows [24]:

3.3.1. Step 1: Initialization Define the design domain and loading conditions. Generate the regular scaled boundary finite element mesh and store the subdomain and nodal information. Define the MIST parameters including the volume fraction of the solid material f solid , the initial element weighting factors xe , the penalization factor b, the initial move limit and the filter radius.

Fig. 13. Dynamic move limit history.

3.3.2. Step 2: The mixed scaled boundary finite element analysis Carry out the mixed SBFEM analysis and calculate the mutual strain energy density. In the proposed technique, the material properties for each subdomain are calculated based on the updated element weighting factor xe and using the material interpolation scheme in Section 3.1. Then in the mixed SBFEM analysis, the element coefficient matrices in Eq. (17) and the solution of the strain and stress in Eqs. (21) and (22) are calculated based on the updated material properties. The objective function is calculated based on the stress and strain and stored.

31

C. Li, L. Tong / Computers and Structures 165 (2016) 24–33

Fig. 15. Pressure contour of the final fluid-rubber mount structures: (a) volume fraction of solid material 45% and (b) volume fraction of solid material 60%.

p=-1

p=-1

p=0 p=-1

p=-1

p=-1

p=0

p=0

p=0

p=-1

p=-1

(a)

(b)

Fig. 16. The pressure fields of the fluid region: (a) volume fraction of solid material 45% and (b) volume fraction of solid material 60%.

3.3.3. Step 3: Generation of the / function and filtering the / surface Construct the surface of the / function by using Eq. (31). In the present study, the mutual strain energy density are first evaluated at the center of each subdomain and then the nodal values are calculated by the Lagrange interpolation method [25]. Finally the / function is constructed. As described in Luo and Tong [25], the 2nd order interpolation functions are used to find the node values over the 3 adjacent subdomains. Mesh dependency problem is common in the topology optimization [22,24]. In order to overcome this issue, the / function is filtered. The values of / are modified based on a linear weighted average within a spatial neighborhood. 3.3.4. Step 4: Find the iso-surface level t and update the element weighting factor xe During the iterative process of topology optimization, an isovalue surface is found to intersect the response function, which forms the structural boundary. The iso-surface threshold value t is obtained by an iterative bisection method when meeting the volume constraint. The area of solid material is summed as PNe i¼1 xeik Ai where Ne is the total number of elements; Ai is the area of element i; k is the current iteration number and xei is the weighting factor of element i. The iso-surface threshold value t is updated as

Ne X t min ðkþ1Þ ¼ t k if xeik Ai > f solid ; t max ðkþ1Þ ¼ t max ðkÞ i¼1

ð32aÞ

or

if

Ne X tmin ðkþ1Þ ¼ t min ðkÞ xeik Ai < f solid ; t max ðkþ1Þ ¼ t k i¼1

ð32bÞ

The iso-surface level t is calculated as

tk ¼

t max ðkÞ þ t min ðkÞ 2

ð33Þ

PNe    xeik Ai The iteration terminates when  f i¼1 Atot  1 < f with a small tolersolid

ance f and total area of elements Atot . The elements with / above the iso-surface at all nodes are the solid material, and those with / below the iso-surface at all nodes are the void [24]. The weighting factors of the elements across the iso-surface are the fraction of the projected area above the isosurface. The updating procedure for element weighting factors is as follows:

xeik

8 þ A > < Aik

þ  A  if  Aik  xei0  6 m i þ  ¼ þ  > : xei0 þ sign Aik  xei0 m if Aik  xei0  > m Ai Ai i

ð34Þ

where xei0 is the element weighting factor before using any bisections. m is the move limit of the weighting factors, which can be used to restrict large variations and stabilize the optimization solution.

32

C. Li, L. Tong / Computers and Structures 165 (2016) 24–33

3.3.5. Step 5: Convergence criteria The optimization solution is converged when one or more convergence parameters are less than the given tolerances such as the change in the weighting factors less than a certain percentage. The convergence parameters are set depending on the problem and the relatively large iterative number is normally chosen to avoid the iteration without ending. The iteration terminates if the convergence criteria are met. Otherwise repeat Steps 2–5. 4. Design examples In order to validate the proposed theoretical framework, four widely studied examples reported in the literature are performed. The compressible and incompressible material designs are presented by considering different Poisson’s ratios. The significant features of the incompressible designs are obtained by comparing with the classical compressible designs. 4.1. Bridge-like structure subjected to one point load This example presents a bridge-like structure subjected to a downward single-point load at the middle point of bottom side. The solid volume fraction is 35% of the total design domain. Fig. 2 shows the geometric configuration and boundary conditions. The same problem was studied in Bruggi and Venini [16]. Plane strain boundary conditions are considered. The difference of the final topologies between compressible and incompressible medias is investigated by considering different Poisson’s ratios. The compressible material property is chosen as Young’s modulus E ¼ 1 and Poisson’s ratio m ¼ 0:25, resulting in K solid ¼ 0:8 and Gsolid ¼ 0:4. The incompressible material property is E ¼ 1 and m  0:5 with K solid ¼ 100 and Gsolid ¼ 0:333. To avoid possible numerical singularity, small values of the void K void ¼ Gvoid ¼ 0:001 are used. The penalization factor b in Eq. (25) equals to 3. A constant move limit of 0.2 is used to stabilize the solution. Due to symmetry of the structure, only half of the design domain is modeled. In order to investigate the effectiveness of the filter, different SBFE meshes are used including 20  20, 40  40 and 80  80 SBFEM subdomains. Correspondingly, the density filter size is chosen as 1.5, 2.0 and 2.5 for three types of mesh. The objective function history versus iteration and the topological boundary contour at iterations 5, 10, 40 and 60 are plotted in Fig. 3 for compressible material using different meshes. The corresponding final compliance is 45.2, 42.7 and 42.1. It is shown the SBFE mesh with 40  40 SBFEM subdomains leads to convergent compliance and topology. Thus the SBFE mesh with 40  40 SBFEM subdomains is used in the subsequent study. The / functions and iso-surfaces for compressible and incompressible design at iterations 5, 10, 40 and 60 are plotted in Figs. 4 and 5, respectively. The area with the / function larger than the iso-threshold value is solid and that with the / function smaller than the iso-threshold value is void. The objective function history versus iteration and the topological boundary contour at iterations 5, 10, 40 and 60 are plotted in Fig. 6 for incompressible design, respectively. The final compliance for incompressible material is C ¼ 35:4. The present result for incompressible material is around 5.8% less than the reference value reported in Bruggi and Venini [16]. Compared to the result for a compressible material, the result for an incompressible material has a higher upper arch, which indicates the structure takes advantage of material incompressibility. 4.2. Bridge-like structure subjected to two point loads In the second example, the same geometry and displacement boundary conditions as the first example are considered. Two point

loads are applied on the middle points of the top and bottom surfaces as shown in Fig. 7. Again, plane strain boundary conditions are considered. The same mesh as in the first example is used. The compressible and incompressible material properties are the same as in Section 4.1. The solid volume fraction is 35% of the total design domain. Fig. 8(a) and (b) shows the comparison between compressible and incompressible material design. Unlike in the first example, the present topology for compressible material is quite different from the one of incompressible material. Fig. 8 shows the objective function history and the topological boundary contour at iterations 5, 10, 40 and 60. The final compliance for incompressible material ðC ¼ 96:3Þ is about 20% less than the one of compressible material ðC ¼ 115:6Þ. The compliance reported in Bruggi and Venini [16] for incompressible design was 100.5. The present result is about 4% less than the reference result. 4.3. Cantilever beam Fig. 9 shows the design domain of a cantilever beam subjected to a symmetric load at the top right corner and lower right corner. The left side of the design domain is fixed. Plane strain conditions are considered. 80  80 SBFEM subdomains are used to model the whole domain. The same problem was studied in Sigmund [29] for compressible material with E ¼ 1 and m ¼ 0:3 (K solid ¼ 0:9615, Gsolid ¼ 0:3846). In the present study, the incompressible material with E ¼ 1 and m  0:5 (K solid ¼ 100, Gsolid ¼ 0:333) is also investigated to show the difference between compressible and incompressible material design. The volume fraction of the solid material is 40%. For the compressible design, the constant move limit of 0.2 and the filter size value of 2.5 are used for iterations. For the incompressible design, the dynamic move limit given in Vasista and Tong [22] is used. The dynamic move limit is used when the objective function starts to oscillate. By using a small value of dynamic move limit, oscillation can be avoided. However, this may lead to slow convergence of objective function. Dynamic move limit is used with the initial value of 0.05 for incompressible material design. It decreases to 0.025 after 90 iterations and increases to 0.05 after 120 iterations. In order to get the converged results, incompressible material design is obtained with a smaller dynamic move limit but with larger iteration number than those of the compressible material design. Furthermore, the dynamic filter size with initial value 3 is applied as shown in Fig. 10. The penalization factor b is 3 in this study. The topological boundary contour at iterations 5, 20 and 40 for compressible design and those at iterations 5, 30, 90 and 150 for incompressible material design are compared and plotted in Fig. 11. The optimized design for the compressible material is a bridge-like shape with six chambers that are symmetric about the middle line. There are two relatively large chambers in the middle region and four small chambers on the lateral sides. However, the incompressible design becomes a five chambers’ topology with a large chamber in the central area and two small chambers on each side. The final compliance for compressible material in Sigmund [29] is C ¼ 17:6. The present result of compressible material design is C ¼ 17:3 which agrees well with the reference result. The present result of incompressible material design is C ¼ 15:4 which is about 10% less than the one of compressible material. 4.4. Fluid-rubber mount structure This example presents a design domain shown in Fig. 12. The exterior thin layer of the structure is fixed to be rubber. The design domain will be filled with rubber and fluid. The uniform pressure with the value of 1 is applied on the top surface. The same problem was studied by Jang and Kim [19]. This problem is analyzed for various volume fraction of rubber within the design domain including

C. Li, L. Tong / Computers and Structures 165 (2016) 24–33

45% or 60%, respectively. The structure is subjected to a distributed pressure on the top surface and the bottom surface is fixed. Again the plane strain assumption is considered. The whole structure is discretized by 80  80 SBFEM subdomains. The interpolation scheme in Eq. (27) is used to represent the rubber or fluid phase. The penalization factor b in Eq. (25) is 1. The dynamic move limit is used to stabilize the solution and plotted in Fig. 13. A density filter size of 2.5 is used. The bulk modulus of the rubber is K solid ¼ 100 and the shear modulus of the rubber is Gsolid ¼ 0:333, whereas the bulk modulus of fluid is K fluid ¼ 100 and the shear modulus is Gfluid ¼ 0:001. Fig. 14(a) and (b) shows the objective function histories with the rubber volume constraint of 45% and 60%, respectively. The final compliance of 45% rubber volume constraint is C ¼ 1:56 and the one with the rubber volume constraint of 60% is C ¼ 1:05. The boundary contour at iterations 1, 10, 50 and 100 with the rubber volume constraint of 45% is shown in Fig. 14(a). Fig. 14(b) shows the boundary contour of the rubber mount at iterations 1, 10, 30 and 60 with the rubber volume constraint of 60%. The final topology of 45% rubber volume constraint consists of six chambers and it is consistent with the one shown in Fig. 10b of Jang and Kim [19]. When the rubber volume fraction increases to 60%, only five chambers exist in the final topology and it shows slightly difference with the one in Fig. 10c of Jang and Kim [19]. The pressure field of the final structures is calculated using Eq. (24). The contours of the pressure fields are plotted in Fig. 15. The pressure fields of the fluid region in the final design structures are shown in Fig. 16 for the two cases. 5. Conclusions In this paper, the semi-analytical mixed SBFEM is extended for topological layout design optimization problems involving incompressible material using the recent developed MIST method. The mixed SBFEM is locking free and passes the inf–sup conditions without requiring any mesh-dependent and stabilization parameters, which enables the topology optimization of incompressible materials possible. The mixed scaled boundary finite element governing equation in displacement is derived and it is solved analytically. This leads to the analytical solutions of stress and strain. Topology optimization method, i.e. MIST, is used to solve a compliant mechanism problem. In MIST, the present response function of the mutual strain energy density is formulated by the stress and strain, which is effective to update the design variables during the iterative process. Numerical studies computed by MIST are presented for the optimal design of compressible and incompressible materials. The present results agree well with the reference results. The numerical results indicate that the MIST topology optimization algorithm based on the mixed SBFEM can yield optimal topologies of both compressible and incompressible materials. The final topologies of incompressible design show significant difference with those of compressible design. Acknowledgment The authors are grateful to the financial support of the Australian Research Council through a Discovery Project Grant (DP140104408).

33

References [1] Herrmann LR. Elasticity equations for incompressible and nearly incompressible materials by a variational theorem. AIAA J 1965;3 (10):1896–900. [2] Malkus DS, Hughes TJ. Mixed finite element methods-reduced and selective integration techniques: a unification of concepts. Comput Meth Appl Mech Eng 1978;15(1):63–81. [3] Hughes TJ. Generalization of selective integration procedures to anisotropic and nonlinear media. Int J Numer Meth Eng 1980;15(9):1413–8. [4] Wilson EL, Taylor RL, Doherty WP, Ghaboussi J. Incompatible displacement models. In: Fenves SJ, Perrone N, Robinson AR, Schnobrich WC, editors. Numerical and computer methods in structural mechanics. New York: Academic Press; 1973. p. 43–57. [5] Simo JC, Rifai MS. A class of mixed assumed strain methods and the method of incompatible modes. Int J Numer Meth Eng 1990;29(8):1595–638. [6] Kasper EP, Taylor RL. A mixed-enhanced strain method: Part I: Geometrically linear problems. Comput Struct 2000;75(3):237–50. [7] Sussman T, Bathe KJ. A finite element formulation for nonlinear incompressible elastic and inelastic analysis. Comput Struct 1987;26(1):357–409. [8] Simo JC, Armero F, Taylor RL. Improved versions of assumed enhanced strain tri-linear elements for 3D finite deformation problems. Comput Meth Appl Mech Eng 1993;110(3):359–86. [9] Fortin M, Brezzi F. Mixed and hybrid finite element methods. Springer; 1991. [10] Bathe KJ. The inf–sup condition and its evaluation for mixed finite element methods. Comput Struct 2001;79(2):243–52. [11] Klaas O, Maniatty A, Shephard MS. A stabilized mixed finite element method for finite elasticity: formulation for linear displacement and pressure interpolation. Comput Meth Appl Mech Eng 1999;180(1):65–79. [12] Lamichhane BP. A stabilized mixed finite element method based on gbiorthogonal systems for nearly incompressible elasticity. Comput Struct 2014;140:48–54. [13] Song C, Wolf J. The scaled boundary finite-element method–alias consistent infinitesimal finite-element cell method–for elastodynamics. Comput Meth Appl Mech Eng 1997;147(3–4):329–55. [14] Li C, Tong L. A mixed SBFEM for stress singularities in nearly incompressible multi-materials. Comput Struct 2015;157:19–30. [15] Sigmund O, Clausen PM. Topology optimization using a mixed formulation: an alternative way to solve pressure load problems. Comput Meth Appl Mech Eng 2007;196(13):1874–89. [16] Bruggi M, Venini P. Topology optimization of incompressible media using mixed finite elements. Comput Meth Appl Mech Eng 2007;196(33):3151–64. [17] Johnson C, Mercier B. Some equilibrium finite element methods for twodimensional elasticity problems. Numer Math 1978;30(1):103–16. [18] Bruggi M, Cinquini C. An alternative truly-mixed formulation to solve pressure load problems in topology optimization. Comput Meth Appl Mech Eng 2009;198(17):1500–12. [19] Jang GW, Kim YY. Topology optimization with displacement-based nonconforming finite elements for incompressible materials. J Mech Sci Technol 2009;23(2):442–51. [20] Panganiban H, Jang GW, Chung TJ. Topology optimization of pressure-actuated compliant mechanisms. Finite Elem Anal Des 2010;46(3):238–46. [21] Jang GW, Panganiban H, Chung TJ. P1-nonconforming quadrilateral finite element for topology optimization. Int J Numer Meth Eng 2010;84 (6):685–707. [22] Vasista S, Tong L. Design and testing of pressurized cellular planar morphing structures. AIAA J 2012;50(6):1328–38. [23] Vasista S, Tong L. Topology-optimized design and testing of a pressure-driven morphing-aerofoil trailing-edge structure. AIAA J 2013;51(8):1898–907. [24] Vasista S, Tong L. Topology optimisation via the moving iso-surface threshold method: implementation and application. Aeronaut J 2014;118 (1201):315–42. [25] Luo Q, Tong L. Design and testing for shape control of piezoelectric structures using topology optimization. Eng Struct 2015;97:90–104. [26] Zienkiewicz OC, Taylor RL. The finite element method, vol. 3. London: McGraw-Hill; 1977. [27] Song C. A matrix function solution for the scaled boundary finite-element equation in statics. Comput Meth Appl Mech Eng 2004;193(23–26):2325–56. [28] Tong L, Lin J. Structural topology optimization with implicit design variable– optimality and algorithm. Finite Elem Anal Des 2011;47(8):922–32. [29] Sigmund O. A 99 line topology optimization code written in Matlab. Struct Multidisc Optimiz 2001;21(2):120–7.