Level set-based topology optimisation of a compliant mechanism design using mathematical programming

We propose a structural optimisation method, based on the level set method and using mathematical programming such as the method of moving asymptotes (MMA), which we apply to the design of compliant mechanisms. A compliant mechanism is a monolithic joint-free mechanism designed to be flexible to obtain a specified motion. In the design of compliant mechanisms, several requirements such as the direction of the deformation and stress concentrations must be considered to obtain the specified mechanical function. Topology optimisation, the most flexible type of structural optimisation, has been successfully used as a design optimisation method for compliant mechanisms, but the utility of topology optimisation results is often spoiled by a plethora of impractical designs such as structures containing grayscale areas. Level set-based topology optimisation methods are immune to the problem of grayscales since the boundaries of the optimal configuration are implicitly represented using the level set function. The proposed method updates the level set function using mathematical programming to facilitate the treatment of constraint functionals. To verify its capability, we apply our method to compliant mechanism design problems that include displacement constraints and stress constraints.


Introduction
Compliant mechanisms are gaining increasing attention as their application in myriad mechanical devices such as MEMS broadens. A compliant mechanism is a monolithic joint-free mechanism designed to be flexible to obtain a specified motion. The major advantages of compliant mechanisms are simplified manufacturing and assembly, reduced cost, lack of mechanical play, silent operation, and freedom from lubrication requirements (Howell, 2001). The first approach to compliant mechanism design was a kinematic synthesis approach in which rigid-body mechanisms were synthesized into compliant mechanisms (e.g., Her and Midha, 1987). This approach, however, is limited to lumped compliant mechanism designs. For the design of fully compliant mechanisms, topology optimisation methods using the continuum synthesis approach are used. In such methods, Sigmund (1997) formulated the objective function as the ratio between input and output forces, called the mechanical advantage. Nishiwaki et al. (1998) presented a structural topol-Correspondence to: M.  ogy optimisation method for compliant mechanisms, where the concept of mutual energy was used in the formulation of flexibility.
Topology optimisation, firstly proposed by Bendsøe and Kikuchi (1988), has been successfully applied to many problems such as minimum mean compliance problems (Suzuki and Kikuchi, 1991), eigen-frequency problems (Diaz and Kikuchi, 1992), electromagnetic problems (Yoo et al., 2000) and so on. The basic concepts of topology optimisation are the extension of the design domain and replacement of the optimisation problem with a material distribution problem using the characteristic function (Murat and Tartar, 1985). Such material distribution problems are known to be illposed, so a relaxation technique is required, such as the homogenization design method or the SIMP method. Topology optimisation methods are the most flexible of optimisation methods since topological changes as well as shape changes are allowed during the optimisation process. However, this advantage is often offset in the optimal configurations by a plethora of impractical designs such as structures containing grayscales or excessive detail, which spoils the utility of the optimal configurations from an engineering standpoint. A type of structural optimisation method using level set boundary expressions has been proposed in which the boundaries of the optimal configuration are implicitly represented using the level set function. A level set-based structural optimisation method was firstly proposed by Sethian and Wiegmann (2000) where the level set function is updated based on the Von Mises stress. Wang et al. (2003) and Allaire et al. (2004) proposed a level set-based structural optimisation method where the level set function is updated using the Hamilton-Jacobi equation, based on the shape sensitivities. Several level set-based structural optimisation methods for the design of compliant mechanisms have been proposed and applied to a multi-material problem (Wang et al., 2005), and a stress minimization problem (Allaire and Jouve, 2008). However, these particular level set-based structural optimisa-tion methods can be categorized as a type of shape optimisation because the introduction of holes is not allowed, but the number of holes can be decreased during optimisation. As a result, the obtained optimal configurations are greatly affected by guesses concerning the initial configuration. To alleviate this problem, Allaire et al. (2005) proposed a level set-based structural optimisation method coupled with the topological gradient method (e.g., Céa et al., 2000).
Level set-based structural topology optimisation methods that do not use a level set function having the property of a signed distance function, resulting in that the introduction of holes is allowed, have also been proposed, such as by Wei and Wang (2009), in which a piecewise constant level set function is used. Their method formulates the objective function as the sum of a primary objective functional and the perimeter of the structure, and a constraint is applied so that the level set function becomes piecewise constant. However, the magnitude of the constraint parameter greatly affects the optimal configuration so that, again, initial configuration settings often determine the utility, or lack thereof, of the obtained optimal configurations. Luo and Tong (2008) proposed a level set-based topology optimisation method incorporating a radial basis function for the design of compliant mechanisms. However, some experience is required when choosing appropriate values for the radial basis function parameters. Yamada et al. (2010) proposed a level set-based topology optimisation method incorporating a fictitious interface energy derived from the phase field concept. In this method, the objective functional is the sum of a primary objective functional and a fictitious interface energy. A piecewise constant level set function is used in this method and the updating scheme uses a reaction-diffusion equation.
The aim of this paper is to extend Yamada's method to the design of compliant mechanisms so that displacement and stress constraints can be easily included. In our proposed method, the level set function is updated using the MMA (Svanberg, 1987) here, to facilitate the treatment of constraint functionals. In Sect. 2, the formulation of the level set-based topology optimisation procedure and optimisation problems are discussed. In Sect. 3, the numerical implementation is discussed, and numerical examples considering displacement and stress constraints are presented in Sect. 4, to confirm the validity and utility of the proposed method.

Level set-based topology optimisation
Here, we briefly discuss the level set-based topology optimisation incorporating a fictitious energy. Topology optimisation is formulated in a fixed design domain D that consists of a domain filled with solid material Ω and a domain filled with void material. Using the characteristic function χ(x) defined as the optimisation problem can be replaced by a material distribution problem. The optimisation problem is formulated as where F is the objective functional, G is a constraint functional and G max is the upper limit of the constraint functional.
The following level set function φ(x) is used to represent the boundaries of the structure, where positive values represent the solid domain, negative values represent the void domain, and zero represents the boundary surfaces.
As a result, the characteristic function χ(x) is replaced by following function, χ φ (x).
The above optimisation problem is an ill-posed problem since the optimal configuration expressed by the level set function is not required to be continuous. To regularize the optimisation problem, Yamada et al. (2010) proposed a regularization technique based on Tikhonov regularization. Thus, the above optimisation problem is replaced with the following: The proposed method updates the level set function φ(x) using mathematical programming, the MMA.

Optimal synthesis of compliant mechanisms
Here, the optimum design of compliant mechanisms is briefly discussed. The main goal of the present optimum design process is to maximize the output displacement in a desired direction. Consider the design domain D where the displacement is fixed at boundary Γ u , an input force t in is applied at boundary Γ in , and a dummy vector t out is introduced at the output port, boundary Γ out , along the desired output direction. Two functions are required for the design of compliant mechanisms. One is to provide sufficient flexibility for deformation along a desired direction specified by a dummy vector t out when an input force is applied. The mutual mean compliance between Γ in and Γ out is used in this paper to formulate the flexibility of the target structure. By maximizing the mutual mean compliance, the output displacement is maximized along the direction of the dummy vector t out . The second requirement is for sufficient stiffness to maintain the integrity of structural shapes when workpiece reaction forces and the input force are applied. Dummy springs are imposed at the input and output ports to represent the input and reaction forces. The optimisation problem under a volume constraint is then formulated as follows.
where u 1 is the displacement when t in is applied at Γ in , V max is the upper limit of the volume constraint, and the other notations are defined as follows.
The sensitivity of the objective functional is simply obtained as follows using the adjoint variable method.
where the adjoint field is defined as follows.

Optimum design problem with mutual mean compliance constraint
The formulation of the optimum design problem is now extended to a problem with a constraint so that the displacement in the direction orthogonal to the desired output direction will be constrained. where l 3 (u 1 ) is the mutual mean compliance when the dummy load t out which is orthogonal to t out is applied, as follows.
The sensitivity of G 2 χ φ is also simply obtained using the adjoint variable method.

Optimum design problem with stress constraint
Here, an optimum design problem with a stress constraint is discussed. Since the utility of compliant mechanisms depends on their structural flexibility, stress concentrations easily occur at thin locations that are subject to repeated flexing. Therefore, the implementation of a stress constraint in the optimisation method can be advantageous for the design of reliable compliant mechanisms that avoid structural failures over their projected lifetime. Several different stress constraint formulations have been studied, such as local stress constraints (e.g. Duysinx and Bendsøe, 1998), global stress constraints (e.g. Martins and Poon, 2005;París et al., 2009) and the block aggregated approach (e.g. París et al., 2010a) which is a hybrid approach combining local and global stress constraints. These formulations are compared in the literature (París et al., 2009;París et al., 2010b). The implementation of local stress constraints uses a straightforward approach, with stress constraints imposed at predefined points such as the centre of finite elements, however this often increases computational demands to the point of intractability. On the other hand, global stress constraints impose a single global constraint that aggregates the effect of all local stress constraints. Although local stress constraints may not be strictly satisfied, the use of a global stress constraint greatly reduces computational demands, and we use a global constraint in this research.
The global stress constraint is formulated as follows. wherẽ where µ is a parameter called the "aggregation parameter". Higher magnitudes of parameter µ impose higher penalties for violated local constraints. σ * in the above equation is defined as follows.
where σ VM represents the Von Mises stress, σ max is an applied stress constraint and ψ e is a parameter called the "stress relaxation coefficient" which is introduced to avoid singularity phenomena and is formulated as follows. where ε is a parameter called the "relaxation factor" which adjusts the magnitude of relaxation and H a (φ) is a Heaviside function to approximate the equilibrium function, which will be described in Sect. 3.3. The sensitivity of the global stress constraints is formulated as follows.
where the adjoint variable is obtained by solving the following equation.
3 Numerical implementation

Optimisation algorithm
The optimisation flowchart is shown in Fig. 1. First, the level set function is initialized. Second, the equilibrium equation is solved using the Finite Element Method (FEM) and the objective functional and constraint functionals are then calculated, also using the FEM. If the objective functional has converged, the optimisation procedure is terminated. If not, the sensitivities of objective and constraint functionals, derived as a continuous expression in the previous section, are computed at Gaussian points of the finite elements. The sensitivities are mapped to the nodes of the finite elements, the level set function is then updated using the MMA, and the process returns to the second step.
Note that the method cited earlier (Yamada et al., 2010), the level set function is updated using a reaction-diffusion equation that is derived based on the Lagrange multiplier method. In the case of multi-constraint problems, the derivation of the Lagrange multiplier used in reaction-diffusion equations becomes complicated, so the proposed method uses the MMA to update the level set function, which facilitates the treatment of constraints. The objective and constraint functionals are approximated using a convex function, and the approximated subproblem is solved at each iteration.

Approximated ∇ 2 φ
In the proposed method, 1 2 τ|∇φ| 2 is used as a regularization term. Therefore, −τ∇ 2 φ must be calculated to compute the sensitivity, and we use the following approximation technique. Note that the variation of the regularization term is formulated in Gurtin (1996). Here, a Dirichlet boundary condition is applied on the non-design boundaries, and Neumann boundary condition is applied on the other boundaries to represent the level set function independently of the outside of the fixed design domain (Yamada et al., 2010). First, we introduce the following time evolutionary equation.
∂D D represents non-design boundaries where the Dirichlet boundary condition is applied. Next, the above equation is discretized in the time domain using the Finite Difference Method, which leads to the following equation.
The above equation is then expressed in weak form as follows.
The above equation can be solved using the FEM. ∇ 2 φ is approximately computed using following equation.

Approximated equilibrium equation
In this paper, the equilibrium equation is approximated using the ersatz material approach, following the literature (e.g. Allaire et al., 2004;Yamada et al., 2010). The equilibrium equation, Eq. (6), is approximated using following equation. where H a (φ) is the Heaviside function defined as follows.
where d is the ratio of the Young's modulus of the void domains to the solid domain and w is the transition width of the Heaviside function.

Numerical examples
In this section, two numerical examples are illustrated to verify the utility and validity of our method. In the following examples, the Young's modulus of the elastic material is 210 GPa, Poisson's ratio is 0.33, and the upper limit of the volume constraint is set to 30% of the design domain. The initial configuration is filled with material in the fixed design domain.
4.1 Compliant mechanism design problem with mutual mean compliance constraint Figure 2 shows the fixed design domain and boundary conditions for model A. The load is applied at the centre of the left edge, with fixed segments at the top of the left edge and extreme left of the bottom edge. A dummy vector t out is applied at the top right of the domain in the horizontal direction. The fixed design domain is descretized using an 80 × 80 mesh of quadrilateral finite elements. The regularization parameter τ is set to 7.0 × 10 −5 . The transition width of the Heaviside function w is set to 0.2 to stabilize the optimisation procedure. Figure 3a shows the optimal configuration without an applied constraint, with Fig. 3c showing its deformed shape. set so that the structure only deforms in the desired horizontal direction, and Fig. 3d shows its deformed shape. As can be seen in the Fig. 3d illustration, the constraint effectively prevents deformation in the orthogonal direction at the output port. The mutual mean compliance represented by dummy vector t out is −2.43 × 10 −9 without a constraint and −1.42 × 10 −10 with a constraint, indicating that displacement in the vertical direction is significantly reduced. The mutual mean compliance represented by dummy vector t out , however, shows little differece in value without and with a constraint, 3.65 × 10 −9 and 3.60 × 10 −9 , respectively. Since the transition width of the Heaviside function w is set to 0.2, the optimal configuration contains grayscale areas to some extent. The grayscale areas in the optimal configuration are removed by setting w = 1.0 × 10 −3 . The mutual mean compliance values represented by dummy vector t out and t out are then −1.42 × 10 −10 and 3.60 × 10 −9 , respectively, which are essentially the same as before removing the grayscale areas.
From the numerical results, we can confirm that the proposed method successfully imposed a displacement constraint for the design of compliant mechanisms, which the existing method, where the level set function is updated using a reaction diffusion equation, cannot easily accomplish. of the left side. Since the design domain is symmetric, only the top half is analyzed in the optimisation process. The applied stress constraint σ max is 5.0×10 3 . The design domain is discretized using an 80 × 80 mesh of quadrilateral finite elements and regularization parameter τ is set to 1.0×10 −4 . The transition width of the Heaviside function w is set to 0.8 and the global stress constraint relaxation factor ε is set to 0.1. Figure 5a shows the optimal configuration without an applied stress constraint and Fig. 5c shows the Von Mises stress distribution of this configuration. Figure 5b shows the optimal configuration with the stress constraint applied, and Fig. 5d shows the corresponding Von Mises stress distribution. The Von Mises stresses at the centre of finite elements are considered. The maximum value of the Von Mises stress is 7.18×10 3 without the stress constraint and 5.83×10 3 with the stress constraint. Although the local stress constraints are not strictly satisfied, the maximum value of the Von Mises stress is reduced and the obtained mutual mean compliance values show little difference, namely, 7.62×10 −10 without the stress constraint and 7.14 × 10 −10 with the stress constraint. Figure 6 shows the density distribution before and after removing grayscale areas of optimal configuration, by setting w = 1.0 × 10 −3 . The maximum value of the Von Mises stress after removing grayscale areas is 5.86 × 10 3 , while the obtained mutual mean compliance is 7.25 × 10 −10 . Although the maximum value of the Von Mises stress is slightly increased after removing grayscale areas, it is sufficiently small compared with the value obtained without an applied stress constraint. Thus, useful optimal configurations can be qualitatively obtained using the proposed method.

Conclusions
This paper proposed a level set-based topology optimisation using mathematical programming. In the proposed method, the level set function is updated using mathematical programming, the MMA, to facilitate the treatment of the constraint functional. This is more difficult with the existing method, in which the level set function is updated using a reaction diffusion equation, where the derivation of the Lagrange multiplier becomes complicated in problems with multiple constraints. A topology optimisation method for compliant mechanisms considering a mutual mean compliance constraint and a stress constraint was presented and optimisation problems were formulated. Although the proposed approach can not explicitly prevent the creation of lumped compliant mechanisms, applying the stress constraint strongly inhibits this, since small areas subject to flex-ure, and notch hinges, tend to be locations where stress is high. The proposed method was applied to compliant mechanism design problems considering a mutual mean compliance constraint and a stress constraint. A global stress constraint is applied but because it does not require that the stress constraint at every point in design domain be satisfied, the optimal configuration does not strictly satisfy all local stress constraints, even though the global stress constraint is satisfied. We confirmed that the maximum stress was effectively reduced in the obtained optimal configuration. Although the optimal configurations contained grayscale areas to some extent, it was confirmed that useful optimal configurations can be qualitatively obtained using the proposed method.