A new level set method for systematic design of hinge-free compliant mechanisms

Comput. Methods Appl. Mech. Engrg. 198 (2008) 318–331

Contents lists available at ScienceDirect

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

A new level set method for systematic design of hinge-free compliant mechanisms Junzhao Luo a, Zhen Luo b, Shikui Chen c, Liyong Tong b,*, Michael Yu Wang c a b c

School of Mechanical Science & Engineering, Huazhong University of Science & Technology, Wuhan, Hubei 430074, China School of Aerospace, Mechanical & Mechatronic Engineering, The University of Sydney, Room 331, Building J11, Sydney, NSW 2006, Australia Department of Mechanical & Automation Engineering, The Chinese University of Hong Kong, Shatin NT, Hong Kong, China

a r t i c l e

i n f o

Article history: Received 11 February 2008 Received in revised form 6 August 2008 Accepted 10 August 2008 Available online 23 August 2008 Keywords: Compliant mechanisms Hinges Topology optimization Level set methods Quadratic energy functionals

a b s t r a c t This paper presents a new level set-based method to realize shape and topology optimization of hingefree compliant mechanisms. A quadratic energy functional used in image processing applications is introduced in the level set method to control the geometric width of structural components in the created mechanism. A semi-implicit scheme with an additive operator splitting (AOS) algorithm is employed to solve the Hamilton–Jacobi partial differential equation (PDE) in the level set method. The design of compliant mechanisms is mathematically represented as a general non-linear programming with a new objective function augmented by the high-order energy term. The structural optimization is thus changed to a numerical process that describes the design as a sequence of motions by updating the implicit boundaries until the optimized structure is achieved under specified constraints. In doing so, it is expected that numerical difficulties such as the Courant–Friedrichs–Lewy (CFL) condition and periodically applied re-initialization procedures in most conventional level set methods can be eliminated. In addition, new holes can be created inside the design domain. The final mechanism configurations consist of strip-like members suitable for generating distributed compliance, and solving the de-facto hinge problem in the design of compliant mechanisms. Two widely studied numerical examples are studied to demonstrate the effectiveness of the proposed method in the context of designing distributed compliant mechanisms. Ó 2008 Elsevier B.V. All rights reserved.

1. Introduction The field of structural shape and topology optimization has recently experienced considerable development reflected by a wide range of engineering applications [1]. In structural topology optimization, a prescribed sum of material is distributed in the design domain subject to specified supports and applied loads via an iterative numerical process until structural performance is extremized. There are many possible applications of the topology optimization method, and the design of compliant mechanisms is one of such potential areas that have attracted much attention in recent years. Compliant mechanisms [2] are a relatively new breed of jointless mechanical devices used to transform displacement, force, or energy from an input port to an output one. A compliant mechanism achieves at least a part of its flexibility from elastic deformation induced by stored strain energy, which in turn leads to a reduction in friction, lubrication, assemblage, noise and vibration. The concept of distributed compliant mechanisms [3] originated from topology optimization using the homogenization method [4], results in the design of various compliant mechanisms with distributed compliance. It is especially suitable for designing * Corresponding author. Tel.: +61 2 93516949; fax: +61 2 93514841. E-mail address: [email protected] (L.Y. Tong). 0045-7825/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2008.08.003

compliant devices in micromechanical systems [5], as some issues associated with fabrication and assemblage can be avoided. Level set method is recently emerged as an alternative approach for performing structural shape and topology optimization. Sethian and Wiegmann [6] are among the first few researchers who introduced the standard level set method [7,8] into structural shape and topology optimization on a fixed Eulerian grid. There are approximately two groups of level set methods that have been developed for solving shape and topology optimization problems. One group is to advance the implicitly represented design boundary by directly solving the Hamilton–Jacobi PDE via explicit schemes [9– 12]. However, it was noted [13,14] that computational difficulties could be encountered in the complicate PDEs solving procedure in the standard level set method. In particular, the final design is dependent on initial guess and can be easily trapped in local minima because of lack of capabilities of generating new holes inside the design domain [15]. Another group is a family of parameterization or equivalent level set methods [13,14,16,17] that does not need to solve the Hamilton–Jacobi PDE of the level set mode using explicit schemes. In the context of designing distributed compliant mechanisms, one additional difficulty is the tendency of forming de-facto hinges in a resulting mechanism [18,19]. The de-facto hinge makes the compliant mechanism analogous to the conventional rigid-body

J. Luo et al. / Comput. Methods Appl. Mech. Engrg. 198 (2008) 318–331

mechanisms, due to the deformation that mainly concentrates on the hinge regions where the mechanism is prone to stress concentration and fatigue breakage. Therefore, there are a number of research attempts to eliminate de-facto hinges in compliant mechanisms. In the framework of material distribution-based methods (i.e., homogenization and SIMP methods), constraints are introduced into the optimization problem to prevent the formulation of hinges from the created mechanism, e.g. the global perimeter constraint method [20], the slope constrained method [21], the local stress constraints scheme [22], the minimum length scale constraint method [23] and the upper-bound constraint approach [24]. Another approach is to employ filtering schemes. Luo et al. [25] devised a density–sensitivity duplicate filtering scheme to control hinges in the final design, and Sigmund [26] presented a morphology-based density filtering scheme for controlling grey scale transitions between solid and void regions so as to generate manufacturable solutions by constraining minimum length scale of structural feature sizes. In addition, Yoon et al. [27] studied a multi-scale topology optimization formulation with translation-invariant embedded wavelet shrinkage. Guest et al. [28] presented a method for imposing a minimum length scale on structural members using nodal variables-based projection approach. Rahmatalla and Swan [29] proposed a hinge-free method via two different sets of artificial spring models. In the context of level set shape and topology optimization methods, for example, Alexandrov and Santosa [30] presented a topology preserving method using a logarithmic barrier penalty term to enforce topological connection of geometrical components. A variational energy approach [31,32] is proposed to control structural geometric width in the context of the standard level set method. Despite many attempts, an efficient and effective scheme to eliminate de-facto hinges is still needed. This paper addresses the de-facto hinge problem of compliant mechanisms by introducing a quadratic energy functional used in image analysis [33] in the level set method. The original objective function is augmented using the high-order energy term that can describe long-range interactions between different sets of contour points without making reference to any particular structural shape, thus leading to hinge-free compliant mechanisms. Furthermore, a semi-implicit scheme with an AOS algorithm [34,35] rather than explicit schemes is used to solve the Hamilton–Jacobi PDE of the level set model so as to avoid numerical difficulties in most standard level set methods.

Fig. 1. Energy changing in the design optimization.

rigid-body mechanism, the input energy will be completely transferred to the output port without energy dissipation. However, what makes a compliant mechanism different from a conventional rigid-body mechanism is that a portion of the input work is stored as elastic energy to enable the ‘‘compliant” movement. For a given input energy, the maximization of output work leads to a minimization of stored elastic energy. As an extreme case, the compliant structure evolves to an analogous rigid-link mechanism (see Fig. 1) when the output work is ideally maximized. In this sense, connections similar to the conventional rigid-body joints, namely the de-facto hinges becomes unavoidable, because only such joints can make a compliant mechanism function as a rigid-link one without energy lose. Thus the de-facto hinges appears unavoidable from the viewpoint of maximizing energy transfer. Hence, this study aims to develop an alternative scheme to design hinge-free complaint mechanisms with modulated or controlled energy transfer. 2.2. Quadratic energy functional As in [31–33], the Euclidean invariant quadratic energy functional is defined as follows:

Eq ðCÞ ¼  2. De-facto hinges and high-order energy functional

¼

2.1. Energy transformation of compliant mechanisms In terms of energy transformation denoted in Fig. 1, the total input energy Ei can be divided into two different parts, namely the output work Eo and the elastic energy Ee stored in the mechanism. The total input energy Ei can be expressed as Ei = FinDin, and the maximal input energy is relevant to the maximal displacement input [Din]max as [Ei]max = Fin[Din]max. The output energy Eo can be defined by

Eo ¼

Z 0

Dout Dgap

F out dx ¼

Z

Dout Dgap

ks x dx ¼ 0

1 ks ðDout  Dgap Þ2 : 2

ð1Þ

And the elastic energy stored in the compliant mechanism is expressed as

1 Ee ¼ Ei  Eo ¼ F in Din  ks ðDout  Dgap Þ2 : 2

ð2Þ

Here, Dgap = f0/ks. The constant ks represents the effective stiffness of the artificial spring mounted at the output port, and f0 refers to a preload applied on the mechanism at the output port. In an ideal

319

¼

Z Z Z Z Z Z

tðpÞ  ~ tðp0 ÞRðj~ dp dp0~ CðpÞ  ~ Cðp0 ÞjÞ ~ tðpÞ ~ tðp0 Þ  Rðj~ CðpÞ  ~ Cðp0 ÞjÞ C p0 j dp0 j~ C p j dpj~ ~ jC p j j~ C p0 j CðpÞ  ~ Cðp0 ÞjÞ: tðpÞ  ~ TðpÞ  ~ Tðp0 ÞRðj~ ds ds0~

ð3Þ

Here, p and p0 are coordinates on the domain of the geometrical curve C that can be expressed in terms of parameter p 2 [0, 1], and ~ tðpÞ ¼ C0ðpÞ is the tangent vector to C at p. ~ TðpÞ and ~ Tðp0 Þ are unit 0 ~ ~ ~ ~ tangent vectors at points CðpÞ and Cðp Þ. CðpÞ and Cðp0 Þ denote the coordinates of two different points on the boundary, respectively, and j~ CðpÞ  ~ Cðp0 Þj denotes the Euclidean distance between point ~ ~ CðpÞ and Cðp0 Þ. ds and ds0 represent the differentials of the arc length. R is a weighting function which was once used to measure the interaction of points within a specified range. In other words, R weight the interactions between different points of the contour according to their distance

8 1; Z if x < dmin  n; > > < h  i xdmin ðxdmin Þ 1 1 ; otherwise; RðxÞ ¼ 2 1  n  p sin p n > > : 0; if x > dmin þ n:

ð4Þ

320

J. Luo et al. / Comput. Methods Appl. Mech. Engrg. 198 (2008) 318–331

Fig. 2. Strip-like shapes controlled by energy functionals with different dmin.

Here, dmin is the radius of the ‘‘hat” function, and n is a half width of the transitional zone. Any points on the contour is related to the parameter dmin of the weighting function R. It should be noted that the function R is always positive. The quadratic energy functional Eq(C) can be regarded as a coherent approach of generating functionals with increased complexity, through which an arbitrary functional can describe arbitrary information about the geometry of a region. It is noted that the quadratic energy is non-local and it can describe interactions or correlations between different points of the contour. The desired shape information is incorporated into the dynamic contour changes via the interactions or correlations, leading to the capability of controlling the geometrical dimension of the final design. Furthermore, non-local forces will be caused by the variational analysis of the quadratic energy term, which in essential indicates the fact that the force at a point of the contour is dependent on both the global boundary shape and its infinitesimal neighborhood through the convolution of the energy along the contour. The mechanism of the quadratic energy functional can be illustrated further in the following discussion. On one hand, a parallel tangent vectors derived from a pair of points interacting with each other will lead to minima of the quadratic term of the energy. The quadratic energy is thus in favor of generating straight boundaries by expanding outwards according to the non-local forces. On the other hand, for a pair of points having anti-parallel tangent vectors, the quadratic component of the energy serves as a softened ‘‘hardcore” [32,33] to make such points mutually repellent, as a result, prevents the points from moving closer than dmin. In this way, the minimization of the quadratic energy tends to lengthen round tips of each region, consequently, likely to generate elongated strip-like structural components with a prescribed width controlled by dmin. In Fig. 2, dmin has a notable influence on the final component width with a mesh size of 30 mm  30 mm. A smaller dmin yields thinner strip-like members with more complex topology while a larger dmin leads to thicker components with relatively simpler topology. Therefore, numerical experience is needed in selecting an appropriate value for dmin to generate desired shape features.

oUðx; tÞ ð5Þ þ vn jrUj ¼ 0; Uðx; 0Þ ¼ U0 ðxÞ; ot where U(x, t) is a Lipschitz continuous scalar function which is used to enable dynamic evolvements of the boundary over a reference domain D  Rd (d = 2 or 3) that includes all admissible shapes of the design domain X(X  D). The motion of the level set in the Hamilton–Jacobi PDE is permitted only along the normal direction of the geometric shape. In the context of diffeomorphism, the change in the embedded geometric shape is only related to the normal component of the velocity field, while its tangential component does not influence the changes of the shape deformation but only its parameterization [8,36]. The normal velocity field is expressed  n = v  ffi (rU/jrUj), v = dx/dt, n = r U/jrUj and as vn =pvffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi jrUj ¼ rU  rU. Hence, moving the design boundary along the normal direction n becomes a question of finding a way to move the scalar function U (x, t) by solving the Hamilton–Jacobi equation on a fixed Eulerian grid. Since a general analytical function for the level set function is usually unknown, time-marching schemes are usually used to solve the Hamilton–Jacobi equation to enable the discrete level set processing. The design domain of compliant mechanisms is illustrated as Fig. 3, where an artificial spring model [37] is applied at the output port to simulate the reaction force from a work-piece. This study uses the mechanical efficiency (ME) to quantify the performance of the compliant mechanism by conducting a well-balanced combination of GA (geometrical advantage) and MA (mechanical advantage) between structural flexibility and mechanical stiffness [14,38]. For simplicity, the present optimization problem only concentrates on linear elasticity but it can be readily extended to geometrically non-linear structures [39–42]. The optimization problem with the augmented objective function is defined as follows:

8 Minimize : > > ðU Þ > > <

Jðu; UÞ ¼ J 0 ðu; UÞ þ kEq ðCÞ 8  > < Ruin  uin 6 0; > > HðUÞdX  V  6 0; > D > Subject to : > : : aðu; v; UÞ ¼ lðv; UÞ; uj

oXd

3. Optimal synthesis of compliant mechanisms In the standard level set method [7,8], the boundary shape of the structure is represented as the zero level set of the level set surface. The Hamilton–Jacobi PDE is defined as

ð6Þ ¼ 0;

8v 2 U;

where J represents the new objective function derived from the original objective functional J0 augmented with the quadratic energy functional Eq(C). The first constraint is introduced to limit the allowable displacement input uin to control the maximal stress level in the resulting mechanism [37], and the second constraint is used

321

J. Luo et al. / Comput. Methods Appl. Mech. Engrg. 198 (2008) 318–331

1 ðuin  uin Þ2 2 K1 Z  Z 2 1 þ k2 HðUÞdX  V  þ HðUÞ  V  : 2 K2 D D

Jðu; UÞ ¼ ½J 0 ðu; UÞ þ kEq ðCÞ þ k1 ðuin  uin Þ þ

ð11Þ

Thus, the shape derivative of the augmented objective function J can be given by

Fig. 3. Compliant mechanisms with artificial spring model.

to generate narrow structures with thin bars and beam-like components by restricting the material usage V*. k is used as a weighting factor to control the energy term in the objective function. The state equation is expressed in terms of the energy bilinear functional a(u, v) and load linear functional l(v).

Z

Dijkl eij ðuÞekl ðvÞHðUÞdX ¼

D

Z

pi vi HðUÞdX Z þ si vi dðUÞjrUjdX;

 oJðu; UÞ o ¼ J 0 ðu; UÞ þ kEq ðCÞ ot ot

o 1 ðuin  uin Þ2 k1 ðuin  uin Þ þ þ ot 2K1 ( Z  Z 2 ) o 1 : HðUÞdX  V  þ HðUÞ  V  k2 þ ot 2 K2 X D ð12Þ Supposing J0 can be expressed in the form of R J0 ðu; UÞ ¼ D f ðuÞHðUÞdX, the shape derivative of J0 is then written as follows:

D

ð7Þ

D

oJ 0 ðu; UÞ ¼ ot

where p and s represent the body forces and boundary tractions, respectively. H(U) is the Heaviside function which serves as a characteristic function to uniformly indicate different parts in the reference domain D and d(U) is the related Dirac function [11]. Using the dummy load method [37], the displacements can be expressed by superposition of two loading cases applied at the input and output ports, respectively. Thus, the objective function J0(u, U) of the mechanical efficiency [14] can be defined as follows:

where

ME ¼ signðGAÞðMA  GAÞ;

aU ðv; wÞ ¼

ð8Þ

where GA is measured by the ratio of displacements Dout and Din at the output and input ports, and MA is the ratio of reaction force Fout and Fin at the output and input ports. The sign(GA) is used to indicate the desired direction of the output displacement. The force Fout = kout  Dout at the output port can be achieved in terms of the spring model, where ks is the stiffness of the artificial spring. GA and MA are defined, respectively, as

GA ¼

Dout u1;out ¼ ; Din u1;in  kout u1;in u2;out þ kout u1;out u2;in

ð9Þ

Z

G0 ðUÞdðUÞjrUjv  n dX;

ð13Þ

D

G0 ðUÞ ¼ f ðuÞ þ pi wi  Dijkl eij ðuÞekl ðwÞ   rU rU þ rðsi wi Þ  si wi : þ r jrUj jrUj

ð14Þ

Here, the adjoint displacement field w can be obtained from the following equation:

Z D

of ðuÞ HðUÞv dX: ou

ð15Þ

Thus, the shape derivative of ME can be expressed as

 oME oGA oMA ; ¼ signðGAÞ MA  þ GA  ot ot ot

ð16Þ

where

!  2 oMA kout ou1;out kout u1;out ou2;out ¼ þ ot 1  kout u2;out ot ot ð1  kout u2;out Þ2

ð17Þ

and

F out kout u1;out ¼ ; MA ¼ F in 1  kout u2;out

ð10Þ

where u1,in, u1,out, u2,in and u2,out are the displacements included in the displacement vectors u1 and u2, respectively. It is noted that u1 and u2 can be obtained by solving two self-adjoint weak forms corresponding to two loading cases f1 and f2, respectively. u1,in and u1,out are displacements at the input and the output ports caused by the input force f1. u2,in and u2,out are the displacements at the input and the output ports caused by applying f2. 4. Shape sensitivity analysis In this section, the concept of the shape derivative [43] is employed as a gradient method to solve the optimization problem in Eq. (6). This paper employs the shape derivative method to measure the sensitivity of boundary perturbations with respect to the pseudo-time. The details of developing the shape derivative can be easily found in the related references [11,12,14]. In this study, the augmented Lagrangian method is applied to convert the original constrained optimization problem into an unconstrained problem as follows:

oGA 1 ou1;out u1;out ðkout u2;out  1Þ ou1;in kout u1;in u1;out ¼ þ þ ot H ot ot H2 H2 ou2;out  ; ot

ð18Þ

where

H ¼ u1;in  kout u1;in u2;out þ kout u1;out u2;in :

ð19Þ

The shape derivative for the objective function J can be further written as follows:

oJ 0 ðu; UÞ oME ¼ ¼ ot ot

Z

G1 ðUÞdðUÞjrUjv  n dX;

ð20Þ

D

where

G1 ðUÞ ¼ f½MAðu; UÞK 1 Dijkl eij ðu1 Þekl ðu1 Þ þ ½MAðu; UÞðK 2 þ K 4 Þ þ GAðu; UÞK 5 Dijkl eij ðu1 Þekl ðu2 Þ þ ½MAðu; UÞK 3 þ GAðu; UÞK 6 Dijkl eij ðu2 Þekl ðu2 Þg;

ð21Þ

where GA and MA can be obtained according to Eqs. (9) and (10), and the six coefficients are defined, respectively, as

322

J. Luo et al. / Comput. Methods Appl. Mech. Engrg. 198 (2008) 318–331

K 1 ¼ ðu1;out  kout u1;out u2;out Þ=ðu1;in  kout u1;in u2;out þ kout u21;out Þ; ð22Þ K 2 ¼ kout u21;out =ðu1;in  kout u1;in u2;out þ kout u21;out Þ;

ð23Þ

K 3 ¼ kout u1;out u1;out =ðu1;in  kout u1;in u2;out þ kout u21;out Þ;

ð24Þ

K 4 ¼ kout u1;in u1;out =ðu1;in  kout u1;in u2;out þ kout u21;out Þ;

ð25Þ

K 5 ¼ kout =ð1  kout u2;out Þ;

ð26Þ

2

K 6 ¼ kout u1;out =ð1  kout u2;out Þ2 :

ð27Þ

Similarly, the shape derivative for the term of Eq(C) can be expressed in terms of [31–33]

) ~ np0 _ Rðj~ rjÞ~ r v  n ds j~ C p0 jdp0 oX oX j~ C p0 j

Z Z ~p0 ds0 v  n ds: _ ~ ¼2 Rðj rjÞ~ rN

o Eq ðCÞ ¼ 2 ot

(I

I

oX

ð28Þ

oX

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~p0 ¼ fxp0 ; y 0 g= x20 þ y20 is the normal Here, r ¼ j~ CðpÞ  ~ Cðp0 Þj and N p p p vector of point p0 . Considering ds = d(U)jrUjdX, the shape derivative of the energy can be given further as

Z oEq ðCÞ ¼ 2 G2 ðUÞdðUÞjrUjv  n dX; ot Z D ~p0 ds0 : _ ~ Rðj rjÞ~ rN G2 ðUÞ ¼

ð29Þ ð30Þ

oX

Thus, the shape derivative of the augmented objective function J can be written as

oJðU; XðtÞÞ oJ 0 ðU; XðtÞÞ oEq ðCÞ ¼ þk ot ot ot Z

GðUÞdðUÞjrUjv  n dX;

¼

The element density and stiffness can be obtained by assuming that they are approximately proportional to the corresponding areafractions of the solid material. 5. Semi-implicit AOS scheme for the level set equation The semi-implicit AOS scheme [35] is used to enable the discrete level set processing. The basic idea behind this method is to split the m-dimensional spatial operator into a set of one-dimensional space discretizations that can be solved using a wellfounded Gaussian elimination algorithm named Thomas Algorithm [35]. The AOS scheme satisfies all discrete scale-space criteria for practical time steps, and it treats all coordinate axes equally and can be easily implemented in arbitrary dimensions. The major features of AOS scheme in topology optimization has been demonstrated in structural mean compliance design [45]. The semi-implicit AOS scheme for updating the level set function is expressed as

Ukþ1 ¼

m 1 X ½I  msbAl ðUk Þ1 ðUk þ sBðUÞÞ; m l¼1

where Bl(Uk) = [I  msAl(Uk)] indicate 1-D diffusion along the xl axes. All these operators are composed of a strictly diagonally dominant tridiagonal linear system. It is noted that the AOS scheme has the same first-order Taylor expansion in s as conventional semi-implicit schemes. Substituting the velocity vN in Eq. (33) into the level set model in Eq. (5) leads to

 dU rU ¼ bjrUjdiv þ jrUjvN : dt jrUj ð31Þ

D

GðUÞ ¼ G1 ðUÞ þ 2kG2 ðUÞ  Z 1 1 þ k1 þ ðuin  uin Þ þ k2 þ HðUÞ  V max : ð32Þ K1 K2 D The multiple Lagrangian multipliers ki(i = 1, 2) related to the volume and the maximal displacement input constraints can be updated using the augmented Lagrangian multiplier method [45]. Ki (i = 1, 2) are set to a small positive value which can be described as Kk+1 = aKk with a < 1. Our numerical tests show that the penalization parameters are stable for structural topology optimization problems if a is chosen between 0.10 and 0.60 for structural mean compliance and compliant mechanism designs. It is noted that the derivative of J agrees with the derivative of J if all the constraints are satisfied. A proper velocity field vN should be selected to ensure the decrease of the objective function. The simplest option is to choose the steepest descent direction [11,12] by letting

vN ¼ GðUÞ:

s

¼ bjrUk j½Al ðUk ÞUkþ1  þ vN jrUk j:

Ukþ1 ¼

m 1 X ½I  msbjrUk jAl ðUk Þ1 ½Uk þ sðvN ðxÞjrUk jÞ m l¼1

G2 ðUÞdðUÞjrUjdX 6 0:

81 2 > 2 k k ; > > hl ðjrUji þjrUjj Þ < aij ¼ 0; P > > 1 2 > : 2 k hl

n2N l ðiÞ

ðjrUji þjrUjkn Þ

j 2 Nl ðiÞ; else;

ð39Þ

; i ¼ j:

The operator A can be spilt into two components in x and y directions (for 2D structures)

ð34Þ

D

In numerical implementation, the discontinuities crossing elements can be approximately evaluated for moving discontinuities using the standard finite element method. The strain energy of an element cut by boundary can be modeled with several methods, e.g. the ‘‘ersatz material” scheme [12], the extended finite element method (X-FEM) [44]. In this study, the ‘‘ersatz material” scheme is adopted to fill the void areas with one relatively weak material.

ð38Þ

with the definition of Al(Uk) = aij(Uk), and aij(Uk) can be given as

ð33Þ

Z

ð37Þ

Hence, with the AOS scheme in Eq. (37), the level set function U can be updated as

Thus, the shape derivative can guarantee a descent direction of the objective function

oJðU; XðtÞÞ ¼ ot

ð36Þ

The formulation in Eq. (36) can be re-written as

Uikþ1  Uki

where

ð35Þ

Fig. 4. Design domain of compliant inverter mechanism.

J. Luo et al. / Comput. Methods Appl. Mech. Engrg. 198 (2008) 318–331

!

A1 ¼

o 1 o ; ox jrUjk ox

!

A2 ¼

o 1 o : oy jrUjk oy

ð40Þ

323

In numerical implementation, the level set function U can be updated as

Fig. 5. Level set contour plots with size control of energy functional term.

Fig. 6. Plots of element stiffness (or density) distribution.

324

J. Luo et al. / Comput. Methods Appl. Mech. Engrg. 198 (2008) 318–331

Fig. 7. Zoom-ins of element-wise stiffness in the hinge regions.

Fig. 8. Contours of local energy density for the mechanism.

Fig. 9. Convergence (0–200 iterations) of the inverter mechanism.

J. Luo et al. / Comput. Methods Appl. Mech. Engrg. 198 (2008) 318–331

Fig. 10. Level set contour plots for mesh-convergence.

Fig. 11. Plots of element stiffness for mesh-convergence with close-ups displays.

Fig. 12. Thresholded energy density distribution contours for the mechanism.

325

326

J. Luo et al. / Comput. Methods Appl. Mech. Engrg. 198 (2008) 318–331

1 2 þ svN ðxÞjrUk j:

Ukþ1 ¼ ½ðI  2sbA1 ðUk ÞÞ1 þ ðI  2sbA2 ðUk ÞÞ1 ½Uk ð41Þ

To simplify the above equation, let

N1 ðUÞ ¼ ðI  2sbA1 ðUk ÞÞ1 ½Uk þ svN ðxÞjrUk j; k

1

k

k

N2 ðUÞ ¼ ðI  2sbA2 ðU ÞÞ ½U þ svN ðxÞjrU j:

ð42Þ ð43Þ

Thus, N1(U) and N2(U) can be easily solved with the following equations using the Thomas method as a Gaussian elimination algorithm for tridiagonal systems.

½I  2sbA1 ðUk ÞN1 ðUÞ ¼ Uk þ svN ðxÞjrUk j; k

k

k

½I  2sbA2 ðU ÞN2 ðUÞ ¼ U þ svN ðxÞjrU j:

ð44Þ ð45Þ

In conventional discrete level set methods with explicit schemes, the time-step size must be sufficiently small to ensure numerical stability due to the consideration of the CFL condition. On a fixed Eulerian grid, only the implicit level set function rather than its partial derivatives is continuous crossing meshes because of the noto-

rious polynomial snaking problem related to the polynomial interpolation [16]. Hence a spatially fine mesh is required to appropriately capture the partial derivatives to avoid numerical artifacts contaminating the solution. However, each time-stepping length in explicit schemes is required to be spatially no more than the minimal grid size [8,36]. A large number of iterations are needed to fully advance the front to convergence via small time-step. It was noted that the present semi-implicit scheme can guarantee stability by satisfying discrete scale-space for any practical time-step sizes in solving the level set equation. Hence it is expected that the AOS scheme performs efficiently because of a relatively large time-step size resulted from the complete relaxation of the CFL condition. It is also well known that the conventional level set methods have a tendency of losing the shape feature of the level set surface and thus lead to too flat and/or too steep gradients due to unwanted numerical diffusion caused by local approximation schemes [8]. In general, it is necessary to regularize the original level set surface using re-initialization scheme to ensure a good approximation of the curvature of the free boundary front [36], so as to resurrect the behavior of the level set surface in the neigh-

Fig. 13. Level set contour for design-dependence on output spring stiffness.

Fig. 14. Plots of element stiffness for design-dependence with zoom-ins displays.

J. Luo et al. / Comput. Methods Appl. Mech. Engrg. 198 (2008) 318–331

bor of the front while keeping the zero level set unchanged. But the error accumulation from re-initializations may move the boundary front if the commonly used iterative procedure is used [46,47]. Furthermore, the global re-initializations prevent the nucleation of new holes inside the material domain [15]. However, the present semi-implicit AOS scheme is capable of generating a globally smooth level set surface. This is because the underlying problem is evaluated over the entire level set function and the steep gradients are readily obtained in the whole fixed and extended design domain. The semi-implicit AOS scheme is unlikely to yield a level set function with lost surface features [35,45], and hence it is capable of creating new holes. It should be pointed out that it is not the AOS scheme that can create new holes, but the elimination of global re-initializations together with the evaluation of the velocity field and design sensitivity over the entire domain. 6. Numerical implementation Two benchmark examples are used to demonstrate the features of the present level set method. The artificial material properties are described as: Young’s modulus for solid material is E = 1 GPa while for the weak material is Emin = 0.001 GPa and Poisson’s ratio is v = 0.31. 6.1. Compliant inverter mechanism The design domain is shown in Fig. 4 with a dimension of 120 lm  120 lm. Only the upper half of the structure is considered with 120  60 finite elements due to its symmetry. The elastic modulus for the artificial material is E = 1 GPa and the Poisson ratio is 0.31. An actuation force Fin = 500 lN is applied at the center point of the left side, and an artificial spring with known stiffness k = 0.1 N/mm is mounted at the output position to simulate the resistance from a work-piece. The maximal material usage is 20% and the allowable input displacement is set as 6 lm. The time-step for the AOS scheme is s = 0.5, the control radius of the quadratic energy functional is dmin = 8 lm and the weighting factor for the energy term is selected as k = 0.15. As indicated by plots of the level set contour at the zero level set in Fig. 5, the quadratic energy term can be used to control the geometric size of the compliant mechanism via an appropriately imposed control radius, leading to a compliant mechanism without the de-facto hinges. To analysis the problem further, the actual element stiffness plots behind the level set contours are shown in Fig. 6 together with zoom-ins of the element-wise stiffness in the hinge regions (Fig. 7). We can find that the regions surrounding the grey scale de-facto hinges are now characterized with distributed solid materials. Also, the topologically optimized contours in

327

Fig. 7 show the distributed local energy density for the mechanism subjected to an input load and an output resistance, in particular, around the neighbour of the hinge regions. It should be pointed out that the white areas inside the structure do not indicate zero energy density distribution but a colour scale with a thresholded energy density (see Fig. 8). It is easy to note that the fixed points on the side near the upper and lower corners and also the input force applied position have a relatively higher compliance. Hence, the compliant mechanism has a more reasonable material and energy density distribution, which are in favour of decreasing stress concentration and possible fatigue breakage. Fig. 9 shows that the convergent curves of the objective function and the volume ratio ranging from 0 to 200 iterations. At the initial 25 iterations, the objective function changes in an opposite direction due to the violation of the volume constraint, and after that the objective function increased to 0.8340 using 223 iterations until the shape and topology converge to a (local) minimum. The curve of the volume constraint denotes that the proposed level set method is mass conservative. It should be pointed out that there is no guarantee that the present level set method can lead to a ‘‘global” optimum because of the non-convexity of the original optimization problem. To study mesh-convergence, the design domain is discretized by using 80  40 = 3200, 120  60 = 7200 and 160  80 = 12 800 elements, respectively. The time-step size is chosen as s = 0.45, the weighting factor for the energy term k = 0.15 and the stiffness for the spring is k = 0.1. The level set contours, element stiffness plots and thresholded energy density contours are shown in Figs. 10–12, respectively. The corresponding values of the mechanical efficiency are 0.8272, 0.8361 and 0.8405, and the iteration numbers at convergence are 223, 267 and 295. The numerical results

Fig. 16. Design domain of the compliant gripper mechanism.

Fig. 15. Thresholded contours of local energy density for the mechanism.

328

J. Luo et al. / Comput. Methods Appl. Mech. Engrg. 198 (2008) 318–331

show that almost the same designs can be obtained with different mesh refinements, and zoom-ins images at hinge regions also de-

note that the final designs are hinge-free mechanisms. For the finite element analysis, a fine mesh is more capable of capturing

Fig. 17. Topologies at different design stages.

Fig. 18. Plots of element stiffness (or density) distribution.

J. Luo et al. / Comput. Methods Appl. Mech. Engrg. 198 (2008) 318–331

the spatial partial derivative of the surface crossing meshing grids, thus to avoid numerical artifacts. But in most material density distribution-based methods, a finer finite element mesh will lead to a more topologically complex structure [1]. Thus, some additional numerical schemes are often incorporated to smear out numerical instabilities such as the checkerboards and the mesh-dependence [20,21]. In most conventional level set methods, a spatially fine mesh indicates a large number of iterations to process the evolution of the level set surface due to the CFL condition, and too many iterations reduce numerical stability due to the accumulation of truncation errors. The example shows that the proposed method appears mesh-independent and also avoids numerical degeneration caused by a large number of iterations. The finite element model with 120  60 = 7200 elements is also used to demonstrate the design-dependence on output spring stiffness. The artificial spring is mounted at output port position with stiffness k = 0.1, k = 0.05 and k = 0.01, respectively. The level set

329

contours, the underlying stiffness plots with zoom-ins and the local energy density contours are given in Figs. 13–15, respectively. It is found that both the mechanical efficiencies (from 0.8245 to 0.8270) and the related convergence iterations (from 229 to 286) remain within a reasonable range. For compliant mechanisms with artificial spring models, the design-dependence on output spring stiffness refers to the fact that the final design relies on the stiffness of the artificial spring attached at the output port position, in general, the smaller of the spring stiffness, and the more complex of the final topology. The numerical results show that the structural topologies are slightly changed with a decrease of the artificial spring stiffness. However, in most conventional methods for distributed compliant mechanisms, small output spring stiffness will often result in a compliant mechanism with de-facto hinges [1]. However, as displayed by the element stiffness and related zoom-in plots in Fig. 14, the present optimization formulation can avoid the hinges even if a relatively small spring

Fig. 19. Zoom-in displays of element-wise stiffness in the hinge regions.

Fig. 20. Contours of local energy density plots.

330

J. Luo et al. / Comput. Methods Appl. Mech. Engrg. 198 (2008) 318–331

Fig. 21. Objective function and volume constraint over 0–250 iterations.

stiffness is employed. Fig. 15 displays the compliant mechanisms produced that are characterized with a relatively uniform energy density distribution represented by a thresholded color contour.

indicates that the design convergences after 297 iterations with a mechanical efficiency of 0.8412, and the volume constraint are mass conservative.

6.2. Compliant gripper mechanism

7. Conclusions

The design domain of the gripper is defined in Fig. 16 with a domain size of 120 lm  120 lm, in which the left-upper and leftlower corners are fixed as the Dirichlet boundary. The lower half of the design domain is discretized with 120  60 finite elements by considering its symmetry, but the topologies of the entire mechanism are displayed. An actuation force Fin = 350 lN is applied at the center point of the left side as a non-homogenous Neumann boundary, and an artificial spring with stiffness k = 0.1 N/mm is attached at the output position to simulate the resistance from the work-piece. The maximal material usage is restricted to 22.5% and the allowable displacement input is limited to 6 lm. The time-step size of the semi-implicit scheme is s = 0.5, the radius in the energy term is dmin = 8 lm and the positive parameter for weighting the energy term in the objective function is given as k = 0.15. The level set contour plots, the element stiffness plots together with zoom-in displays of the hinge regions and the contours of local energy density are depicted in Figs. 17–19, respectively. The level set topologies in Fig. 17 show that the final mechanism can be geometrically described with a smooth boundary and the shape fidelity and topological changes can be performed at the same time, which should be the advantages of the level set-based free boundary representation method. The numerical procedure shows that the present algorithm allows a large time-step size with improved computational efficiency. New holes can be inserted inside the design domain via the boundary merging and splitting processes partially because the global re-initialization procedure is eliminated. Evidently, several numerical difficulties in the conventional level set methods are avoided in this study. Furthermore, the final design in Fig. 17f shows a hinge-free compliant mechanism as a result of introducing the energy functional term. The quadratic energy functional encourages formulating straight boundaries by expanding outwards in terms of the non-local forces and simultaneously preventing the points from moving too closer, and generating elongated strip-like structural components. The underlying plots of element stiffness with the grey scale transition region (Fig. 18) and the related zoom-in displays in hinge areas (Fig. 19) also demonstrate that the present method can generate compliant mechanisms with distributed local energy density (Fig. 20). Fig. 21

This study proposes a new level set method for shape and topology optimization of distributed compliant mechanisms. A quadratic energy functional from image active contour techniques is introduced into the level set model to augment the original objective function, which can lead to hinge-free compliant mechanisms with distributed compliance by controlling structural shape features. The semi-implicit AOS scheme used in image processing is introduced to solve the optimization problem. The AOS scheme is stable for practical time-step sizes and free from the restriction of the CFL condition compared to explicit schemes, because all the terms relevant to stability are discretized in an implicit manner. The numerical examples show that the present method is capable of designing hinge-free compliant mechanisms with controlled strip-like shapes and eliminating numerical difficulties in most conventional level set methods. Acknowledgements This research was funded in part by the Australian Research Council (ARC-DP0666683), the Research Grants Council of Hong Kong (CUHK416507). References [1] M.P. Bendsøe, O. Sigmund, Topology Optimization: Theory, Methods, and Applications, Springer, Berlin, Heidelberg, 2003. [2] L.L. Howell, Compliant Mechanisms, John Wiley & Sons Inc., New York, 2001. [3] G.K. Ananthasuresh, S. Kota, Y. Gianchandani, A methodical approach to the design of compliant micromechanisms, in: Solid-State Sensor and Actuator Workshop, 1994, pp. 189–192. [4] M.P. Bendsøe, N. Kikuchi, Generating optimal topology in structural design using a homogenization method, Comput. Methods Appl. Mech. Engrg. 71 (1988) 197–224. [5] G.K. Ananthasuresh, L.L. Howell, Mechanical design of compliant microsystems – a perspective and prospects, J. Mech. Des. 127 (2005) 736–738. [6] J.A. Sethian, A. Wiegmann, Structural boundary design via level set and immersed interface methods, J. Comput. Phys. 163 (2000) 489–528. [7] S. Osher, J.A. Sethian, Front propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations, J. Comput. Phys. 78 (1988) 12–49. [8] J.A. Sethian, Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer version and material science, Cambridge Monograph on Applied and Computational Mathematics, Cambridge University Press, Cambridge, UK, 1999.

J. Luo et al. / Comput. Methods Appl. Mech. Engrg. 198 (2008) 318–331 [9] S. Osher, F. Santosa, Level-set methods for optimization problem involving geometry and constraints: I. Frequencies of a two-density inhomogeneous drum, J. Comput. Phys. 171 (2001) 272–288. [10] G. Allaire, F. Jouve, A.M. Toader, A level-set method for shape optimization, CR Math. 334 (2002) 1125–1130. [11] M.Y. Wang, X.M. Wang, D.M. Guo, A level set method for structural topology optimization, Comput. Methods Appl. Mech. Engrg. 192 (2003) 227–246. [12] G. Allaire, F. Jouve, A.M. Toader, Structural optimization using sensitivity analysis and a level-set method, J. Comput. Phys. 194 (2004) 363–393. [13] S.Y. Wang, M.Y. Wang, Radial basis functions and level set method for structural topology optimization, Int. J. Numer. Methods Engrg. 65 (2006) 2060–2090. [14] Z. Luo, L.Y. Tong, M.Y. Wang, S.Y. Wang, Shape and topology optimization of compliant mechanisms using a parameterization level set method, J. Comput. Phys. 227 (2007) 680–705. [15] M. Burger, B. Hacker, W. Ring, Incorporating topological derivatives into level set methods, J. Comput. Phys. 194 (2004) 344–362. [16] T. Belytschko, S.P. Xiao, C. Parimi, Topology optimization with implicitly function and regularization, Int. J. Numer. Methods Engrg. 57 (2003) 1177– 1196. [17] H. Haber, A multilever, level-set method for optimizing eigenvalues in shape design problems, J. Comput. 198 (2004) 518–534. [18] J. Diaz, O. Sigmund, Checkerboard patterns in layout optimization, Struct. Multidisciplin. Optim. 10 (1995) 40–45. [19] O. Sigmund, J. Petersson, Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima, Struct. Multidisciplin. Optim. 16 (1998) 68–75. [20] R.B. Haber, M.P. Bendsøe, C. Jog, A new approach to variable topology shape design using a constraint on the perimeter, Struct. Multidisciplin. Optim. 11 (1996) 1–12. [21] J. Petersson, O. Sigmund, Slope constrained topology optimization, Int. J. Numer. Methods Engrg. 41 (1998) 1417–1434. [22] P. Duysinx, M.P. Bendsøe, Topology optimization of continuum structures with local stress constraints, Int. J. Numer. Methods Engrg. 43 (1998) 1453–1478. [23] T.A. Poulsen, A new scheme for imposing a minimum length scale in topology optimization, Int. J. Numer. Methods Engrg. 57 (2003) 741–760. [24] L. Yin, G.K. Ananthasuresh, Design of distributed compliant mechanisms, Mech. Based Des. Struct. Machines 31 (2001) 151–179. [25] Z. Luo, L.P. Chen, J.Z. Yang, Y.Q. Zhang, K. Abdel-Malek, Compliant mechanism design using a multi-objective topology optimization scheme of continuum structures, Struct. Multidisciplin. Optim. 30 (2005) 142–154. [26] O. Sigmund, Morphology-based black and white filters for topology optimization, Struct. Multidisciplin. Optim. 33 (2007) 401–424. [27] G.H. Yoon, Y.Y. Kim, M.P. Bendsøe, O. Sigmund, Hinge-free topology optimization with embedded translation-invariant differentiable wavelet shrinkage, Struct. Multidisciplin. Optim. 27 (2004) 139–150. [28] J.K. Guest, J.H. Prévost, T. Belytschko, Achieving minimum length scale in topology optimization using nodal design variables and projection functions, Int. J. Numer. Methods Engrg. 61 (2004) 238–254.

331

[29] S. Rahmatalla, C.C. Swan, Sparse monolithic compliant mechanisms using continuum structural topology optimization, Int. J. Numer. Methods Engrg. 62 (2005) 1579–1605. [30] O. Alexandrov, F. Santosa, A topology-preserving level set method for shape optimization, J. Comput. Phys. 204 (2005) 121–130. [31] S.K. Chen, PhD Dissertation, Compliant mechanisms with distributed compliance and characteristic stiffness: A level set approach. Department of Mechanical and Automation Engineering, The Chinese University of Hong Kong, Shatin, NT, Hong Kong, China, 2007. [32] S.K. Chen, M.Y. Wang, A.Q. Liu, Shape feature control in structural topology optimization, Computer Aided Design (2008), doi:10.1016/j.cad.2008.07.004. [33] M. Rochery, I. Jermyn, J. Zerubia, Higher order active contours, Int. J. Comput. Vision 69 (2006) 27–42. [34] T. Lu, T. Neittaanmaki, X.C. Tai, A parallel splitting up method and its application to Navier–Stokes equations, Appl. Math. Lett. 4 (1991) 25–29. [35] J. Weickert, B.M. Romeny, M. Viergever, Efficient and reliable schemes for nonlinear diffusion filtering, IEEE Trans. Image Process. 7 (1998) 398–410. [36] S. Osher, R.P. Fedkiw, Level Set Methods and Dynamic Implicit Surface, Springer-Verlag, New York, 2002. [37] O. Sigmund, On the design of compliant mechanisms using topology optimization, Mech. Struct. Machines 25 (1997) 493–524. [38] M.Y. Wang, S.K. Chen, X.M. Wang, Y.L. Mei, Design of multi-material compliant mechanisms using level set methods, J. Mech. Des. 127 (2005) 941–956. [39] O. Sigmund, Design of multiphysics actuator using topology optimization – part I: one material structure, Comput. Methods Appl. Mech. Engrg. 190 (2001) 6577–6604. [40] C.B.W. Pedersen, T. Buhl, O. Sigmund, Topology synthesis of largedisplacement compliant mechanisms, Int. J. Numer. Methods Engrg. 50 (2001) 2683–2705. [41] T.E. Bruns, D.A. Tortorelli, Topology optimization of nonlinear elastic structures and compliant mechanisms, Comput. Methods Appl. Mech. Engrg. 190 (2001) 3443–3459. [42] Z. Luo, L.Y. Tong, A level set method for shape and topology optimization of large-displacement compliant mechanisms, Int. J. Numer. Methods Engrg. (2008), doi:10.1002/nme.2352. [43] J. Sokolowski, J.P. Zolesio, Introduction to Shape Optimization: Shape Sensitivity Analysis, Springer, Berlin, 1992. [44] T. Belytschko, C. Parimi, N. Moës, N. Sukumar, S. Usui, Structured extended finite element methods for solids defined by implicit surfaces, Int. J. Numer. Methods Engrg. 56 (2002) 609–635. [45] J.Z. Luo, Z. Luo, L.P. Chen, L.Y. Tong, M.Y. Wang, A semi-implicit level set method for structural shape and topology optimization, J. Comput. Phys. 227 (2008) 5561–5581. [46] D. Peng, B. Merriman, S. Osher, H. Zhao, M. Kang, A PDE-based fast local level set method, J. Comput. Phys. 155 (1999) 410–438. [47] M. Sussman, E. Fatemi, An efficient interface-preserving level set re-distancing algorithm and its application to interfacial incompressible fluid flow, SIAM J. Sci. Comput. 20 (1999) 1165–1191.