An augmented formulation of distributed compliant mechanism optimization using a level set method

Abstract

Topology optimization has emerged as one of the key approaches to design compliant mechanisms. However, one of the main difficulties is that the resulted compliant mechanisms often have de facto hinges. For this reason, a simple yet efficient formulation for designing hinge-free compliant mechanisms is developed and examined within a level set–based topology optimization framework. First, the conventional objective function is augmented using an output stiffness. Second, the proposed formulation is solved using a level set method for designing some benchmark problems in the literature. It is shown that the proposed augmented objective function can prevent the de facto hinges in the obtained compliant mechanisms. Finally, some concluding remarks and future work are put forward.

Introduction

A compliant mechanism is regarded as a mechanism that gains its mobility from the deflection of its flexible members.^1,2 Over the past decades, compliant mechanisms have been extensively studied. Several approaches have been developed for the design of compliant mechanisms. Generally, these approaches can be categorized into two types. The first one is the kinematics synthesis approach. ¹ In this approach, the compliant mechanism is derived from a known rigid-body mechanism.^1,3 Although the method has been successfully used in designing compliant mechanisms for precision applications, it requires a good deal of designers’ intuition and involvement.

The second approach is derived from the structural topology optimization approach.^2,4–7 As one of the most challenging tasks in the optimization design, topology optimization has been deeply explored and applied to a variety of design problems, for example, the minimum mean compliance problem, ⁸ the vehicle component design problems, ⁹ and the top-down structural assembly synthesis problem. ¹⁰ During the last decades, several methods have been developed, such as the ground structure method,^11–13 the solid isotropic material with penalization (SIMP) method,^8,14 and the level set method.^15–17

At the very beginning, the level set method was developed for numerically tracking fronts and free boundaries.^18,19 After Sethian and Wiegmann ²⁰ first introduced it to design structural boundaries, Osher and Santosa ²¹ extended the method by introducing the shape sensitivity analysis into the framework. In level set–based topology optimization methods, the structural boundary is treated as the design parameter and is implicitly embedded into a scalar function and updated by solving the level set equation. ²² Therefore, topology changes, for example, merging and splitting, can be easily handled. Furthermore, several numerical instabilities that usually occurred in the density-based topology optimization approaches, such as checkerboard patterns and gray scales, can be eliminated^23–25 as well.

When using continuum topology optimization methods to design compliant mechanisms, one of the significant challenges is their strong tendency to result in de facto hinges. ⁴ In the context of material distribution–based topology optimization methods, several methods have been developed to eliminate the de facto hinges in the design of compliant mechanisms. It is surely possible to redesign the de facto hinge regions as continuous material bridges. ²⁶ However, the obtained mechanisms will deviate from the original mechanism. Alternatively, one may use the procedures derived from the image processing method, ²⁷ such as the filtering method, ²⁸ to eliminate the de facto hinges. These methods can insure the absence of the one-node connected hinges. However, the lumped compliance, which is caused by the de facto hinges, is sometimes inevitable. Researchers also tried to develop new design models, such as Rahmatalla and Swan ⁴ developed a new spring-based method, to prevent de facto hinges. Another method can be regarded as a so-called hybrid discretization method.^29,30 In this method, during each optimization iteration, each design element will be subdivided for finite element analysis. This often leads to a low computational efficiency although the de facto hinges can be eliminated.

For using the level set method, a logarithmic barrier penalty term has been utilized to ensure topological connection of the structural components. ³¹ However, this cannot eliminate the de facto hinges. A possible strategy is to control the geometric width of structural components using a quadratic energy functional, as stated in Luo et al. ³² and Chen et al. ³³ In addition, an intrinsic characteristic stiffness method is developed by Wang and Chen. ³⁴ Deepak et al. ³⁵ indicate that this method can also lead to point flexures when a large objective geometrical advantage (GA) is needed. Zhu and Zhang ⁶ developed two alternative formulations for developing hinge-free compliant mechanisms using level set method. This work has been further extended for designing compliant mechanisms with multiple outputs. ⁷

The fundamental reason for the occurrence of the de facto hinges lies in the mathematical formulation of the design problem. ³⁶ For this reason, in this study, we were trying to propose a new objective function for the design of hinge-free compliant mechanisms and verified its validity by designing several two-dimensional (2D) numerical examples that are widely studied in the literature of compliant mechanism optimization.

Level set method

In level set method, the structural boundary $\partial \Omega$ is implicitly embedded in a scalar function ϕ^18,19 as its zero level set. Therefore

$$\{\begin{array}{ccc}\varphi (\mathbf{x},t)>0\hfill & \mathrm{if}\hfill & \mathbf{x}\in \Omega \hfill \\ \multicolumn{1}{c}{\varphi (\mathbf{x},t)=0}& \mathrm{if}\hfill & \mathbf{x}\in \partial \Omega \hfill \\ \multicolumn{1}{c}{\varphi (\mathbf{x},t)<0}& \mathrm{if}\hfill & \mathbf{x}\in D\backslash \Omega \hfill \end{array}$$

(1)

where D is the reference domain to contain all permissible shapes of the design domain $\Omega$ , $\partial \Omega$ is the interface of the structure, and $D\backslash \Omega$ represents the void area.

The optimization is achieved by solving the following Hamilton–Jacobi equation

$$\frac{\partial \varphi }{\partial t}+V_n|\nabla \varphi |=0$$

(2)

where t is the time and velocity $V_n$ determines the motion of the interface.^15,16

Optimal synthesis of compliant mechanisms and the challenge of de facto hinges

Conventional mathematical formulation of the optimization problem

Synthesis of compliant mechanisms has been formulated in many different ways under two main groups. The first of which is to establish the objective function by maximizing a mechanical measurement, such as the mechanical advantage (MA), ² the GA,^6,37 and the mechanical efficiency (ME). ³² The second of which is formulated by considering both flexibility and compliance to meet the function and strength requirements. ¹¹ A comparison study of these formulations can be found in Deepak et al. ³⁵

For designing compliant mechanisms with single input–output behavior, the design domain can be illustrated in Figure 1 where ${\Gamma }_d$ indicates the Dirichlet boundary. An input force $f_{\mathrm{in}}$ is applied at the input port i. $u_{\mathrm{in}}$ and $u_{\mathrm{out}}$ indicate the displacements occurred at the input port i and the output port o due to $f_{\mathrm{in}}$ , respectively. A spring with stiffness $k_{\mathrm{out}}$ is attached to the output port to imitate the reaction force from the workpiece by $f_{\mathrm{out}}=k_{\mathrm{out}}{u}_{\mathrm{out}}$ .

Figure 1.

The design domain of topology optimization of compliant mechanisms.

In this study, in selecting the stiffness of $k_{\mathrm{out}}$ , the bounding spring value method ⁴ is employed. The value of $k_{\mathrm{out}}$ is set to be ${10}^{-4}{k}_b$ where $k_b$ is computed by applying a unit load at the input port of the reference domain D when D is fully occupied by the structural material and the output port is restrained, see Rahmatalla and Swan. ⁴

Here, we choose GA to quantify the performance of the compliant mechanisms. Incorporating with the level set method, a conventional formulation for topology optimization of the compliant mechanisms can be formulated as follows ³⁷

$$\begin{array}{cc}min:\hfill & -\mathrm{GA}(u,\varphi )=-\frac{u_{\mathrm{out}}}{u_{\mathrm{in}}}\hfill \\ \multicolumn{1}{c}{\mathrm{s}.\mathrm{t}.:}& u_{\mathrm{in}}\le u_{\mathrm{in}}^{\max }\hfill \\ \multicolumn{1}{c}{}& \mathrm{Vol}=\underset{D}{\int }H(\varphi )d\Omega \le \mathrm{Vo}{l}^{\max }\hfill \\ \multicolumn{1}{c}{}& a(u,v,\varphi )=l(v,\varphi )\forall v\in U\hfill \end{array}$$

(3)

where $u_{\mathrm{in}}$ is constrained by an upper limit $u_{\mathrm{in}}^{\max }$ for indirectly controlling the maximum stress level. ² Vol denotes the total material usage and is constrained by an upper limit $\mathrm{Vo}{l}^{\max }$ . v denotes the virtual displacement fields in space U. Note that a minus is used in equation (3) since maximization of GA is equivalent to minimization of $-\mathrm{GA}$ . The energy bilinear functional $a(u,v,\varphi )$ and the load linear functional $l(v,\varphi )$ are, respectively, expressed as follows

$$a(u,v,\varphi )=\underset{D}{\int }{E}_{\mathrm{ijkl}}{\epsilon }_{\mathrm{ij}}(u){\epsilon }_{\mathrm{kl}}(u)H(\varphi )d\Omega$$

(4)

$$l(v,\varphi )=\underset{D}{\int }\mathrm{fv}\delta (\varphi )|\nabla \varphi |d\Omega$$

(5)

where $E_{\mathrm{ijkl}}$ and ${\epsilon }_{\mathrm{ij}}$ are the material property tensor and strain tensor, respectively. Since the design condition considered in this article does not concern the body force, only the boundary traction f is considered in the above equations.

Using the dummy load method, ² the displacements $u_{\mathrm{out}}$ and $u_{\mathrm{in}}$ can be expressed by superposition of two loading cases applied at the input and output ports. Therefore

$$\mathrm{GA}=\frac{u_{1,o}}{u_{1,i}-k_{\mathrm{out}}{u}_{1,i}{u}_{2,o}+k_{\mathrm{out}}{u}_{1,o}{u}_{2,i}}$$

(6)

where $u_{1,i},u_{1,o},u_{2,i}$ and $u_{2,o}$ are the displacements included in the displacement field $u_1$ and $u_2$ , respectively. $u_1$ is the displacement field which is obtained by applied a unit load $f_1$ at the input port of the design domain. $u_2$ is the displacement field which is obtained by applied a unit load $f_2$ at the output port of the design domain. $u_{1,i}$ and $u_{1,o}$ are the displacements at the input and the output ports, respectively, caused by $f_1$ . $u_{2,i}$ and $u_{2,o}$ are the displacements at the input and the output ports, respectively, caused by $f_2$ . For more details, please refer to Sigmund ² and Chen. ³⁷

$H(\varphi )$ is the Heaviside function defined as follows

$$H(\varphi )=\{\begin{array}{ccc}\begin{array}{c}1\hfill \\ \multicolumn{1}{c}{0}\end{array}\hfill & \begin{array}{c}\mathrm{if}\hfill \\ \multicolumn{1}{c}{\mathrm{if}}\end{array}\hfill & \begin{array}{c}\varphi \ge 0\hfill \\ \multicolumn{1}{c}{\varphi <0}\end{array}\hfill \end{array}$$

(7)

and $\delta (\varphi )$ is the one-dimensional delta function defined as follows

$$\delta (\varphi )=\frac{\mathrm{dH}(\varphi )}{d\varphi }$$

(8)

De facto hinge problem

A hinged compliant inverter which is obtained using equation (3) is shown in Figure 2. These kinds of hinges are not needed since they make the obtained compliant mechanisms very difficult to fabricate, especially in the micro-scale. ⁴ Since the flexibility of the mechanism is only provided in localized areas (hinge areas), the stress in the hinge areas would approach very high and the mechanism would break.

Figure 2.

A compliant displacement inverter mechanism suffers from the de facto hinges, which are marked with dashed line circle.

The reason for de facto hinges lies behind the objective formulations. ³⁶ In fact, for designing compliant mechanism, many developed formulations specify two main purposes, that is, maximizing the elastic deformation at the output port, meanwhile minimizing the overall compliance. This makes the true optimum of the problem a rigid-body linkage with revolute joints. The reason is that it can generate the largest output motion and has the minimum strain energy. From this point of view, the de facto hinges are inevitable. There are some other formulations that try to avoid these two main purposes to avoid the de facto hinges, such as the characteristic stiffness formulation. ³⁷ A comparative review of those formulations can be found in Deepak et al. ³⁵

The output stiffness and a new formulation

For topological synthesis of a compliant mechanism, although the optimum mechanism could not be known in advance, all the three regions (input region, output region, and fixed region) must be connected to one another in order to form a meaningful structure. ⁵ Here, we proposed an output compliance that can be built in the conventional formulation (equation (3)). The idea of designing the hinge-free compliant mechanisms by augmenting the conventional optimization problem with additional energy functionals is not new. For example, Luo et al. ³² proposed an augmented objective function using a quadratic energy functional to control the geometric width of the mechanism. However, it is difficult to implement this formulation to other topology optimization methods, such as the SIMP method or evolutionary structural optimization method. ³⁸

As shown in Figure 3, the output stiffness $E_{\mathrm{out}}$ is determined based on the case that an external unit force is only applied at the output port o while keeping the input port i as a free boundary (unfixed). Therefore, $E_{\mathrm{out}}$ can be illustrated as follows

$$E_{\mathrm{out}}=\underset{\partial D}{\int }{f}_2{u}_{2}d\Gamma$$

(9)

and this makes the proposed formulation very easy to use since no extra finite element analysis needs to be addressed.

Figure 3.

Schematic for determining $E_{\mathrm{out}}$ .

The reason for introducing the energy function to prevent de facto hinges can be stated as follows. As shown in Figure 4(a), for a compliant mechanism that suffers de facto hinges, when the mechanism is loaded with a force at the output port, the corresponding displacement at the output port will approach very high since the surrounding materials of the de facto hinges undergo essentially rigid-body rotations. This is also fairly intuitive as shown in Figure 2. However, when the mechanism is completely free of de facto hinges as is shown in Figure 4(b), when the mechanism is loaded with a force at the output port, the corresponding displacement at the output port will become smaller, that is

$$u_h>u_{\mathrm{hf}}$$

(10)

Figure 4.

Schematic of a (a) hinged and (b) hinge-free compliant mechanism.

This means a hinged compliant mechanism corresponds to a large displacement while a hinge-free compliant mechanism corresponds to a small one. Conversely, by minimizing the displacement at the output port due to the unit load (which is equivalent to minimize the output stiffness $E_{\mathrm{out}}$ ), one can obtain a compliant mechanism which is free of the de facto hinges.

Based on the above analysis, the objective function can be set by minimizing $-\mathrm{GA}$ meanwhile minimizing $E_{\mathrm{out}}$ . Incorporating the level set model, a new mathematical formulation of the optimization problem can be rewritten as follows

$$\begin{array}{cc}min:\hfill & -\mathrm{GA}(u,\varphi )+\alpha E_{\mathrm{out}}(u,\varphi )\hfill \\ \multicolumn{1}{c}{\mathrm{s}.\mathrm{t}.:}& u_{\mathrm{in}}(\varphi )\le u_{\mathrm{in}}^{\max }\hfill \\ \multicolumn{1}{c}{}& \underset{D}{\int }H(\varphi )d\Omega \le \mathrm{Vo}{l}^{\max }\hfill \\ \multicolumn{1}{c}{}& a(u,v,\varphi )=l(v,\varphi )\forall u\in U\hfill \end{array}$$

(11)

where $\alpha$ is the weighting factor of $E_{\mathrm{out}}$ . The new objective function derived from the original function (equation (3)) is augmented with the energy function $E_{\mathrm{out}}$ . The minimization of $E_{out}$ is to make the mechanism structure stiffer, while the maximization of GA is to make the mechanism more flexible.

Shape sensitivity analysis

In order to obtain the sensitivity of the objective function with respect to the boundary perturbations, the shape derivative method^39,40 is employed. The optimization problem represented in equation (11) is reformulated using Lagrange’s method of undetermined multipliers as follows

$$\begin{array}{cc}\multicolumn{1}{c}{J=}& -\mathrm{GA}+\alpha E_{\mathrm{out}}+{\lambda }_{u_{\mathrm{in}}}(u_{\mathrm{in}}-u_{\mathrm{in}}^{max}+{\theta }_{u_{\mathrm{in}}}^2)\hfill \\ \multicolumn{1}{c}{}& +{\lambda }_{\mathrm{Vol}}(\mathrm{Vol}-\mathrm{Vo}{l}^{max}+{\theta }_{\mathrm{Vol}}^2)\hfill \end{array}$$

(12)

where ${\theta }_{u_{\mathrm{in}}}$ , ${\theta }_{\mathrm{Vol}}$ are slack variables to convert the inequality constraint into the equality one and ${\lambda }_{u_{\mathrm{in}}}$ , ${\lambda }_{\mathrm{Vol}}$ are Lagrange multipliers.

Applying the Kuhn–Tucker conditions of J leads to

$$\begin{array}{cc}\multicolumn{1}{c}{〈\frac{\partial J}{\partial \varphi },\phi 〉=}& -〈\frac{\partial \mathrm{GA}}{\partial \varphi },\phi 〉+\alpha 〈\frac{\partial E_{\mathrm{out}}}{\partial \varphi },\phi 〉\hfill \\ \multicolumn{1}{c}{}& +{\lambda }_{u_{\mathrm{in}}}〈\frac{\partial u_{\mathrm{in}}}{\partial \varphi },\phi 〉+{\lambda }_{\mathrm{Vol}}〈\frac{\partial \mathrm{Vol}}{\partial \varphi },\phi 〉=0\hfill \end{array}$$

(13)

$$\frac{\partial J}{\partial {\lambda }_{u_{\mathrm{in}}}}=u_{\mathrm{in}}-u_{\mathrm{in}}^{max}+{\theta }_{u_{\mathrm{in}}}^2=0$$

(14)

$$\frac{\partial J}{\partial {\lambda }_{\mathrm{Vol}}}=\mathrm{Vol}-\mathrm{Vo}{l}^{max}+{\theta }_{\mathrm{Vol}}^2=0$$

(15)

$$\frac{\partial J}{\partial {\theta }_{u_{\mathrm{in}}}}=2{\lambda }_{u_{\mathrm{in}}}{\theta }_{u_{\mathrm{in}}}=0$$

(16)

$$\frac{\partial J}{\partial {\theta }_{\mathrm{Vol}}}=2{\lambda }_{\mathrm{Vol}}{\theta }_{\mathrm{Vol}}=0$$

(17)

$${\lambda }_{u_{\mathrm{in}}}\ge 0,{\lambda }_{\mathrm{Vol}}\ge 0$$

(18)

where $〈\partial J/\partial \varphi ,\phi 〉$ indicates the Fréchet derivative of J with respect to $\varphi$ in the direction of $\phi$ .

Equation (13) can be expressed as

$$〈\frac{\partial J}{\partial \varphi },\phi 〉=\underset{D}{\int }\left(\begin{array}{c}-S_{\mathrm{GA}}+\alpha S_{E_{\mathrm{out}}}\hfill \\ \multicolumn{1}{c}{+{\zeta }_{u_{\mathrm{in}}}{S}_{u_{\mathrm{in}}}+{\zeta }_{\mathrm{Vol}}{S}_{\mathrm{Vol}}}\end{array}\right)\delta (\varphi )\phi d\Omega$$

(19)

where $S_{\mathrm{GA}}$ , $S_{E_{\mathrm{out}}}$ , $S_{u_{\mathrm{in}}}$ , and $S_{\mathrm{Vol}}$ can be derived from $〈\partial \mathrm{GA}/\partial \varphi ,\phi 〉$ , $〈\partial E_{\mathrm{out}}/\partial \varphi ,\phi 〉$ , $〈\partial u_{\mathrm{in}}/\partial \varphi ,\phi 〉$ , and $〈\partial \mathrm{Vol}/\partial \varphi ,\phi 〉$ , respectively. And ${\zeta }_{u_{\mathrm{in}}}$ , ${\zeta }_{\mathrm{Vol}}$ can be determined from the mentioned Kuhn–Tucker optimally condition as follows

$${\zeta }_{u_{\mathrm{in}}}=\{\begin{array}{ccc}{\lambda }_{u_{\mathrm{in}}}\hfill & \mathrm{if}\hfill & u_{\mathrm{in}}>u_{\mathrm{in}}^{max}\hfill \\ \multicolumn{1}{c}{0}& \mathrm{if}\hfill & u_{\mathrm{in}}\le u_{\mathrm{in}}^{max}\hfill \end{array}$$

(20)

$${\zeta }_{\mathrm{Vol}}=\{\begin{array}{ccc}{\lambda }_{\mathrm{Vol}}\hfill & \mathrm{if}\hfill & \mathrm{Vol}>\mathrm{Vo}{l}^{max}\hfill \\ \multicolumn{1}{c}{0}& \mathrm{if}\hfill & \mathrm{Vol}\le \mathrm{Vo}{l}^{max}\hfill \end{array}$$

(21)

In order to calculate the shape sensitivity of J, we need to calculate the shape sensitivities of the individual functions, that is, $〈\partial \mathrm{GA}/\partial \varphi ,\phi 〉$ , $〈\partial E_{\mathrm{out}}/\partial \varphi ,\phi 〉$ , $〈\partial u_{\mathrm{in}}/\partial \varphi ,\phi 〉$ , and $〈\partial \mathrm{Vol}/\partial \varphi ,\phi 〉$ .

Taking Fréchet derivative of the GA with respect to $\varphi$ in the direction of $\phi$ leads to

$$\begin{array}{c}〈\frac{\partial \mathrm{GA}}{\partial \varphi },\phi 〉=\frac{\partial \mathrm{GA}}{\partial u_{1,i}}〈\frac{\partial u_{1,i}}{\partial \varphi },\phi 〉+\frac{\partial \mathrm{GA}}{\partial u_{1,o}}〈\frac{\partial u_{1,o}}{\partial \varphi },\phi 〉\hfill \\ +\frac{\partial \mathrm{GA}}{\partial u_{2,i}}〈\frac{\partial u_{2,i}}{\partial \varphi },\phi 〉+\frac{\partial \mathrm{GA}}{\partial u_{2,o}}〈\frac{\partial u_{2,o}}{\partial \varphi },\phi 〉\hfill \end{array}$$

(22)

Furthermore, taking Fréchet derivative of the $u_{\mathrm{in}}$ with respect to $\varphi$ in the direction of $\phi$ leads to

$$\begin{array}{cc}\multicolumn{1}{c}{〈\frac{\partial u_{\mathrm{in}}}{\partial \varphi },\phi 〉=}& \frac{\partial u_{\mathrm{in}}}{\partial u_{1,i}}〈\frac{\partial u_{1,i}}{\partial \varphi },\phi 〉+\frac{\partial u_{\mathrm{in}}}{\partial u_{1,o}}〈\frac{\partial u_{1,o}}{\partial \varphi },\phi 〉\hfill \\ \multicolumn{1}{c}{}& +\frac{\partial u_{\mathrm{in}}}{\partial u_{2,i}}〈\frac{\partial u_{2,i}}{\partial \varphi },\phi 〉+\frac{\partial u_{\mathrm{in}}}{\partial u_{2,o}}〈\frac{\partial u_{2,o}}{\partial \varphi },\phi 〉\hfill \end{array}$$

(23)

Since $\partial \mathrm{GA}/\partial u_{1,i}$ , $\partial \mathrm{GA}/\partial u_{1,o}$ , $\partial \mathrm{GA}/\partial u_{2,i}$ , $\partial \mathrm{GA}/\partial u_{2,o}$ , $\partial u_{\mathrm{in}}/\partial u_{1,i}$ , $\partial u_{\mathrm{in}}/\partial u_{1,o}$ , $\partial u_{\mathrm{in}}/\partial u_{2,i}$ , and $\partial u_{\mathrm{in}}/\partial u_{2,o}$ can be directly obtained from equation (6), only $〈\partial u_{1,i}/\partial \varphi ,\phi 〉$ , $〈\partial u_{1,o}/\partial \varphi ,\phi 〉$ , $〈\partial u_{2,i}/\partial \varphi ,\phi 〉$ , and $〈\partial u_{2,o}/\partial \varphi ,\phi 〉$ need to be solved.

The Fréchet derivative of $u_{1,i}$ with respect to $\varphi$ in the direction of $\phi$ can be written as

$$\begin{array}{cc}〈\frac{\partial u_{1,i}}{\partial \varphi },\phi 〉=\hfill & \underset{D}{\int }\left(\begin{array}{c}2[\kappa f_1{u}_1+\frac{\partial (f_1{u}_1)}{\partial n}]\hfill \\ \multicolumn{1}{c}{-E_{\mathrm{ijkl}}{\epsilon }_{\mathrm{ij}}(u_1){\epsilon }_{\mathrm{kl}}(u_1)}\end{array}\right)\delta (\varphi )\phi d\Omega \hfill \\ \multicolumn{1}{c}{}& +2{\int }_{\Gamma }\frac{\delta (\varphi )}{|\nabla \varphi |}\frac{\partial \varphi }{\partial n}{f}_1{u}_1\phi d\Gamma \hfill \end{array}$$

(24)

where $\kappa$ is the mean curvature of the structural boundary defined by $\kappa =\mathrm{div}·n$ , where n is the normal vector of the structural boundary. For design cases where the external force is applied at a point, equation (24) is reduced to

$$〈\frac{\partial u_{1,i}}{\partial \varphi },\phi 〉=-\underset{D}{\int }{E}_{\mathrm{ijkl}}{\epsilon }_{\mathrm{ij}}(u_1){\epsilon }_{\mathrm{kl}}(u_1)\delta (\varphi )\phi d\Omega$$

(25)

Similarly, the Fréchet derivatives of $u_{1,o}$ , $u_{2,i}$ , $E_{\mathrm{out}}$ , and $u_{2,o}$ can be written as

$$\begin{array}{c}\multicolumn{1}{c}{〈\frac{\partial u_{1,o}}{\partial \varphi },\phi 〉=〈\frac{\partial u_{2,i}}{\partial \varphi },\phi 〉}\\ \multicolumn{1}{c}{=\underset{D}{\int }{E}_{\mathrm{ijkl}}{\epsilon }_{\mathrm{ij}}(u_1){\epsilon }_{\mathrm{kl}}(u_2)\delta (\varphi )\phi d\Omega }\end{array}$$

(26)

$$\begin{array}{c}\multicolumn{1}{c}{〈\frac{\partial E_{\mathrm{out}}}{\partial \varphi },\phi 〉=〈\frac{\partial u_{2,o}}{\partial \varphi },\phi 〉}\\ \multicolumn{1}{c}{=-\underset{D}{\int }{E}_{\mathrm{ijkl}}{\epsilon }_{\mathrm{ij}}(u_2){\epsilon }_{\mathrm{kl}}(u_2)\delta (\varphi )\phi d\Omega }\end{array}$$

(27)

The derivative of allowable material usage Vol with respect to $\varphi$ in the direction $\phi$ can be written as follows

$$〈\frac{\partial \mathrm{Vol}}{\partial \varphi },\phi 〉=\underset{D}{\int }\frac{\mathrm{dH}(\varphi )}{d\varphi }\phi d\Omega =\underset{D}{\int }\delta (\varphi )\phi d\Omega$$

(28)

For more details, please refer to previous studies.^16,39–41

Numerical implementation

For the implementation of the proposed method, a number of numerical issues need to be addressed here. The second-order accurate essentially non-oscillatory (ENO2) ¹⁸ is employed for solving equation (2).

In order to ensure stability, the following re-initialization equation ⁴² is used

$$\frac{\partial \varphi }{\partial t}+S({\varphi }_0)(|\nabla \varphi |-1)=0$$

(29)

where $S({\varphi }_0)$ is a sign function taken as 1 in $\Omega$ , $-1$ in $D\backslash \Omega$ , and 0 on the interface. In numerical implementation, $S(\varphi )$ is approximated by

$$S({\varphi }_0)=\frac{{\varphi }_0}{\sqrt{{\varphi }_0^2+{(\Delta x)}^2}}$$

(30)

where $\Delta x$ is set to the mesh size.

In order to avoid the numerical difficulties, we adopted the method proposed in Wang et al. ¹⁵ in which the Heaviside function is approximated using equation

$$H(\varphi )=\{\begin{array}{ccc}ϵ\hfill & \mathrm{if}\hfill & \varphi <-\Delta \hfill \\ \multicolumn{1}{c}{(\frac{3}{4}-\frac{3ϵ}{4})(\frac{\varphi }{\Delta }-\frac{{\varphi }^3}{3{\Delta }^3})+\frac{1+ϵ}{2}}& \mathrm{if}\hfill & -\Delta \le \varphi \le \Delta \hfill \\ \multicolumn{1}{c}{1}& \mathrm{if}\hfill & \varphi >\Delta \hfill \end{array}$$

(31)

where $ϵ$ is small value to ensure the numerical stiffness nonsingular. In our case, $ϵ$ is set to be ${10}^{-3}$ . $\Delta$ is set to be $0.75\Delta x$ . Using equations (31) and (8), $\delta (\varphi )$ can be obtained directly as

$$\delta (\varphi )=\{\begin{array}{ccc}0\hfill & \mathrm{if}\hfill & \mathrm{others}\hfill \\ \multicolumn{1}{c}{(\frac{3}{4}-\frac{3ϵ}{4})(\frac{1}{\Delta }-\frac{{\varphi }^2}{{\Delta }^3})}& \mathrm{if}\hfill & -\Delta \le \varphi \le \Delta \hfill \end{array}$$

(32)

Numerical examples

In this section, in order to demonstrate the validity of the proposed formulation, we present several examples. The used artificial material properties are as follows: Young’s modulus for solid material is $E=1$ and Poisson’s ratio is $\upsilon =0.3$ . The void area is assumed with a Young’s modulus $E=0.001$ and the same Poisson’s ratio $\upsilon =0.3$ .

Displacement inverter

The designing of the displacement inverter is first considered. The design domain is shown in Figure 5. The length and height of the design domain have the same size. The left top corner and the left bottom corner are fixed. A single horizontal load $F=1$ is applied at the center point of the left side of the design domain. An inverse displacement is expected at the center point of the right side of the design domain. Due to symmetry, only half the design domain is taken into consideration. $80\times 40$ finite elements are used for discretization.

Figure 5.

The design domain of the displacement inverter topology optimization problem.

Effect of $\alpha$ on the optimum topology

This section is focused on examining the effectiveness of $E_{\mathrm{out}}$ on preventing the de facto hinges. The effect of $\alpha$ on the optimal topology is examined. A volume ratio $\mathrm{Vo}{l}^{\max }=0.3$ is considered. For all studied cases, $u_{\mathrm{in}}^{\max }$ is set to 50.

Five cases are studied, in which $\alpha$ is set to 0.001, 0.002, 0.003, 0.004, and 0.006. The initial design as well as the corresponding final designs are shown in Figure 6. For all studied five cases, the created compliant mechanisms do not suffer de facto hinges. When large $\alpha$ is used, the elastic hinges can also be eliminated and compliant mechanism that only contain strip-like members is obtained which is in favor of generating distributed compliance such as case of $\alpha =0.006$ .

Figure 6.

Topology optimization of the displacement inverter using equation (11) with different $\alpha$ : (a) the initial configuration, (b) $\alpha =0.001$ , (c) $\alpha =0.002$ , (d) $\alpha =0.003$ , (e) $\alpha =0.004$ , and (f) $\alpha =0.006$ .

The value of $\alpha$ has direct impact on the GA of the created mechanism, see Figure 7. In fact, one can find that a larger $\alpha$ surely will lead to a small GA. A proper selection of $\alpha$ depends on applications, and no universal optimal setting may exist. However, to reduce the difficulty of selection, one may rephrase the objection function in (equation (11)) using the normalization method so that the value of $\alpha$ lies between 0 and 1.

Figure 7.

The effect of $\alpha$ on the geometric advantages and the iterations of the obtained mechanisms.

In fact, equation (11) can prevent de facto hinges not only in the final topology but also during the optimization process. This can be seen from some intermediate designs of the displacement inverter with $\alpha =0.002$ in Figure 8. It should be pointed out that the present method cannot generate new holes freely inside the design domain. Therefore, for using the proposed method, a good initial guess of the topology (with a certain number of holes) is needed for obtaining a reasonable solution. The convergence histories of using the proposed method with the case of $\alpha =0.002$ are shown in Figure 9. It needs more than 190 iterations before convergence. Future works will be focused on reducing the computational effort without affecting the optimization outcomes.

Figure 8.

The intermediate designs of the displacement problem with $\alpha =0.002$ : (a) step 1, (b) step 25, (c) step 50, (d) step 100, (e) step 150, and (f) step 190.

Figure 9.

The convergence histories of the displacement inverter problem with $\alpha =0.002$ .

Push gripper

The design of the push gripper has been widely studied previously, and the design domain and boundary conditions are shown in Figure 10. The goal of the design is to achieve a mechanism that when a horizontal force is applied at the input port of the mechanism, the opposing output ports move vertically and therefore it is capable of gripping a workpiece. The design domain is discretized by $80\times 80$ finite elements for the elastic analysis. The void area (gap size) is set to be $30\times 30$ finite elements. Due to symmetry, only half of the design domain is taken into consideration. The weighting factor $\alpha$ is set to 0.004. The maximum material usage $\mathrm{Vo}{l}^{\max }$ is set to $30\%$ and the $u_{\mathrm{in}}^{\max }$ is set to 60.

Figure 10.

The design domain of the push gripper topology optimization problem.

The final designs obtained using the proposed formulation are shown in Figure 11. The corresponding level set surface plots are shown in Figure 12. Note that only half of the design domain is plotted in Figure 12 due to the symmetry. One can confirm that there are no de facto hinges occurring in the created mechanisms and this can confirm the capability of the proposed formulation (equation (11)) for designing hinge-free compliant mechanisms.

Figure 11.

The final designs of the push gripper mechanism.

Figure 12.

Level set surfaces of the optimal topologies of the push gripper mechanism (only half of the design domain is plotted).

Conclusion

A simple yet efficient formulation for topology optimization of hinge-free compliant mechanisms is presented by taking into consideration an output stiffness. The level set method is used for modeling the optimization problem. The proposed method is examined by topology optimization of two benchmark compliant mechanisms, that is, the displacement inverter and the push gripper. It is shown that the augmented objective formulation can prevent de facto hinges. Although only the level set method is considered, the implementation of the proposed formulation to other topology optimization methods should be straightforward. The present method can be used as an alternative method to partially control the de facto hinges of the linear elastic mechanisms. Our future research will investigate the validity of the presented method for designing both 2D and three-dimensional (3D) large-displacement compliant mechanisms.