Design of Large-Displacement Compliant Mechanisms by Topology Optimization Incorporating Modified Additive Hyperelasticity Technique

This paper is focused on the topology design of compliant mechanisms undergoing large displacement (over 20% of the structural dimension). Based on the artificial springmodel and the geometrically nonlinear finite element analysis, the optimization problem is formulated so as tomaximize the output displacement under a givenmaterial volume constraint. Amodified additive hyperelasticity technique is proposed to circumvent numerical instabilities that occurred in the low-density or intermediate-density elements during the optimization process. Compared to the previous method, the modified technique is very effective and can provide more accurate response analysis for the large-displacement compliant mechanism. The whole optimization process is carried out by the gradient-based mathematical programming method. Numerical examples of a force-inverting mechanism and a microgripping mechanism are presented. The obtained optimal solutions verify the applicability of the proposed numerical techniques and show the necessity of considering large displacement in the design problem.


Introduction
With the advantages of fast fabrication and easy miniaturization, compliance mechanisms have been widely used in modern drivetrain systems, aerospace structures, electronic equipment, bioengineering, biomedical devices, and many other applications.The optimal layout of compliant mechanisms can be typically found by using the topology optimization methods.Such a design problem can be viewed as a trade-off between the flexibility to deliver required motions and the rigidity to overcome external resistance.
Early attention to the topology optimization of compliant mechanisms was focused on the mathematical modeling of the design problem.Among others, Ananthasuresh et al. [1] and Saxena and Ananthasuresh [2] used the combination of output displacement and strain energy as the objective function.In this formulation, an actuation force was applied at the input port and a spring was added to the output port.In order to make the optimal mechanisms-like structure adapted to the stiffness of the workpiece, Sigmund [3] presented an optimization formulation to maximize the mechanical advantage (MA) of compliant mechanisms.Hetrick and Kota [4] proposed an energy efficiency criterion, which is the ratio between the net energy transferred at the output and the input energy.In addition, an artificial spring formulation [5] has also been introduced to address this problem, in which the actuator's stiffness was modeled by an input spring and only the output displacement was taken as the objective function.Recently, a comparative study by Deepak et al. [6] revealed that all above formulations yield almost similar topologies, though they show different convergence ability.Based on these optimization formulations and the linear finite element analysis, topology optimization of smalldisplacement compliant mechanisms has been successfully realized by using some familiar approaches including the solid isotropic material with penalization (SIMP) approach [7], the evolutionary structural optimization (ESO) method [8,9], the level set-based topology optimization method [10,11], and some nodal variable-based methods [12,13].
In many circumstances, compliant mechanisms undergo large displacement during operation.To obtain a meaningful optimal topology, it is necessary to accurately describe the nonlinear large-displacement structural behavior of compliant mechanisms in the optimization process.Therefore, many works have been published in the field of topology optimization of compliant mechanisms by incorporating the geometrically nonlinear finite element analysis.Based on the SIMP approach and the total Lagrangian finite element formulation, Pedersen et al. [14] proposed a synthesis tool for the design of large-displacement and path-generating compliant mechanisms.Their results showed that the practical output displacement obtained by nonlinear topology optimization would be up to 2.5 times the output displacement obtained by the linear method.Sigmund [15] further extended the nonlinear method to the topology optimization of thermoelectrical micro actuators.Bruns and Tortorelli [16] investigated the topology optimization problem of compliant mechanisms considering geometrical and material nonlinearities.They suggested an element removal and reintroduction strategy to relieve the instability associated with mesh distortion in lowdensity elements [17].Maute and Frangopol [18] designed microelectromechanical systems (MEMS) by geometrically nonlinear topology optimization accounting for stochastic uncertainties.Luo and Tong [19] developed a parameterization level set method for shape and topology optimization of compliant mechanisms involving large deformation.In these studies, the motion of compliant mechanisms is typically limited by about 10% of the length scale of the design domain.
As the literature survey reveals, topology optimization with geometrically nonlinear analysis provides an attractive way for the design of compliant mechanisms.However, few studies address the topology optimization of compliant mechanisms undergoing very large displacement (e.g., over 20% of the structural dimension), partly because of the extreme distortion of mesh and serious "element instability phenomenon" exhibited in the optimization process.In this case, the tangent stiffness matrix of distorted elements will become negative definite, which leads to serous nonconvergence of the nonlinear finite element analysis as well as the topology optimization process.
In this paper, we aim to develop an effective method of the topology optimization for large-displacement compliant mechanisms.The artificial spring formulation associated with the material interpolation scheme is adopted in the optimization model.The additive hyperelasticity technique reported in our previous work [20] is extended and improved to not only eliminate the element instability phenomenon that often appears in problems with large deformation, but also provide more accurate response analysis.After obtaining the sensitivity information, the topology optimization problem is solved by using the Method of Moving Asymptotes (MMA) [21].Finally, two numerical examples are presented to validate the proposed technique and to investigate the effects of large displacement on the optimal topology of compliant mechanisms.

Void Solid
Figure 1: Design problem of compliant mechanisms.

Topology Optimization Formulation of Compliant Mechanisms
We consider a general design problem of compliant mechanisms for maximizing specified output displacement or output force under a given input force.As mentioned in the Introduction, several optimization formulations have been proposed for such a problem.Here, the artificial spring model proposed by Bendsøe and Sigmund [5] is employed.To this end, two artificial springs are attached to the design domain.
As shown in Figure 1, one spring is applied at the output port to simulate the resistance from a workpiece, and another is put at the loading point to exert a certain control on the net input work.By taking the material relative densities as design variables, the topology optimization of compliant mechanisms is formulated as max   out () where  out () denotes the resulting output displacement,  = [ 1 ,  2 , . . .,   ]  is the vector of elemental relative densities,  is the total finite element number of the discretized design domain, and u and R are the global displacement vector and the global residual force vector, respectively.Note that R depends nonlinearly on the elemental material properties, which are modeled as a power-law function of the elemental relative densities by the SIMP approach.The equilibrium state of the large-displacement compliant mechanism is represented by R(u(), ) = 0 and should be solved by an iterative procedure, for example, the Newton-Raphson method.  is the volume of the eth element and  * is the upper bound on material volume. min = 0.001 is the lower bound of material relative density.

Modified Additive Hyperelasticity Technique.
During the course of material density-based topology optimization, lowdensity elements are inevitable, especially when a sensitivity filter or density filter technique is applied.For problems involving large deformation, these low-density elements may cause severe instability due to local buckling or mesh distortions.The local instability induced by low-density areas may cause serious convergence issue in the geometrically nonlinear finite element analysis.Indeed, this phenomenon is one of the main difficulties hindering the topology optimization of large-displacement compliant mechanisms.Two simple and intuitive treatments, the "convergence criterion relaxation" [22] and the "element removal" [17], have been proposed to overcome local instability in the minimum-density elements.The former method excludes the nodes surrounded by minimum-density elements from the convergence criterion, while the latter method removes minimum-density elements from the finite element mesh.However, both methods cannot completely overcome this difficulty in the design of compliant mechanisms undergoing a very large deformation due to the fact that local instability may also occur in intermediate-density elements (e.g.,   = 0.01 to 0.3).Other more complex methods, such as the element connectivity parameterization method [23,24] and the Levenberg-Marquardt method [25], have also been introduced to stabilize the optimization process.In addition, Wang et al. [26] proposed a new energy interpolation scheme to alleviate the numerical instability in the low stiffness region.In this study, a modified method based on the additive hyperelasticity technique proposed by Luo et al. [20] is employed.
The basic idea of the additive hyperelasticity technique is to add a specified soft hyperelastic material to the low-density elements that may become unstable.The additive hyperelastic material for the eth element has the following strain energy function according to the Yeoh model [27]: where  1 = tr(C) is the first invariant of the right Cauchy-Green strain tensor C,   1 > 0 and   2 > 0 are material constants of the additive hyperelastic material for the eth element,   is the elemental density, and  is the penalization factor used in the SIMP approach.It is seen that the shearing capacity of the additive hyperelastic material is gradually increased as a result of increasing strain.Thus, the stiffness of low-density element is enhanced under compression to prevent itself from instability.In order to achieve this, two issues need to be further addressed.
The first issue is how to identify the low-density elements that tend to be unstable during the optimization process.A simple criterion is the ratio between the element strain in the kth optimization step and the specified threshold value  * ; that is, where  ,von is the average von Mises strain of the eth element.
From our numerical experiments, we set the threshold value as  * = 1.0 in the compliant mechanism topology optimization problem.If  ()  ≥ 1, this means that the deformation of the eth element is close to the instability point.Otherwise, there is little possibility for the eth element to become unstable in the next iteration step.
The second issue is the choice of parameters   1 and   2 for the additive hyperelastic material.Under small strain, local instability is not an issue and thus initial Young's modulus ( ini = 6  1 ) of the additive hyperelastic material should be small enough.Therefore, we set the minimum value   1 =  min   0 /6, where  0 is Young's modulus of the solid elastic material.Under large strain, the stiffness of the additive hyperelastic material highly depends on the parameter   2 and a large value of   2 can ensure effective elimination of instability.However, it may also lead to a large prediction discrepancy of the actual behavior of the mechanism.A reasonable value of   2 for the additive material in the eth element can be updated by an adaptive strategy, which is expressed by where the superscripts ( + 1) and () denote the iteration steps of optimization.Compared with our previous work [20], which uses a uniform value of additive material parameters for all elements, this study proposed a modified strategy to eliminate the element instability by selecting a more reasonable value of   2 for each element.As will be illustrated in Section 4.1, by using the proposed modified algorithm, the discrepancy between the responses of the remodeled structure and the actual structure for large-displacement compliant mechanisms will be substantially reduced.

Nonlinear Finite Element Analysis
Procedures.We now perform the finite element analysis of the remodeled compliant mechanism undergoing large displacement and rotations.The remodeled mechanism is composed of original linear elastic material and additive hyperelastic material.Here, the additive hyperelastic material is invoked to describe nonlinear material behaviors under large strains.Using linear Hooke's law, we express the second Piola-Kirchhoff stress S ori in the original elastic material as follows: where  is the Green-Lagrange strain and D ori is the constitutive elasticity tensor of the original elastic material.Using the SIMP model, for each element, the elasticity tensor D ori  is a function of the elemental relative density   : where D 0 is the elasticity tensor of the fully solid material and  > 1 is the penalization factor.In this study, we set  = 3 for penalizing intermediate-density values during the optimization process.
For the additive hyperelastic part, the second Piola-Kirchhoff stress S add can be obtained from the derivative of the strain energy function Ψ add ; that is, where  1 ,  2 , and  3 are the three invariants of C.Then, according to (2), the second Piola-Kirchhoff stress in the additive hyperelastic material for each element is computed by Under large displacement, the Green-Lagrange strain is defined by and can be rewritten after discretization as where B 0 is the linear part of the strain-displacement matrix and B  represents the second-order term.The Green-Lagrange strain is work-conjugate to the second Piola-Kirchhoff stress.Therefore, the static equilibrium state of the remodeled mechanism is expressed by where R denotes the residual force, Ω 0 is the undeformed volume of the design domain, K spring is the stiffness matrix obtained by assembling the contributions of two artificial springs, and P is the vector of prescribed input force.The static equilibrium (11) for the compliant mechanism is nonlinear and can be typically solved by using the Newton-Raphson procedure.

Design Sensitivity Analysis.
In this section, the adjoint variable method for sensitivities calculation of the objective function  out () is presented.Direct differentiation of the static equilibrium function (11) Then, the gradient  out / is computed by where R/u = K  is the global tangent stiffness matrix at the equilibrium state.
To calculate (13), we first introduce an adjoint vector , which can be easily obtained as the solution to the linear adjoint equation The partial derivative of the global residual with respect to the eth design variable is Here, the partial derivatives of the second Piola-Kirchhoff stresses S ori and S add with respect to the design variable are obtained from the differentiation of ( 5) and (8).Submitting ( 14) and ( 15) into ( 13), it yields that In the finite element implementation, the 4-node plane stress elements are used.It is well known that the use of loworder finite elements in topology optimization may lead to checkerboard patterns and mesh dependency.In this study, the sensitivity filter technique [28] is employed to circumvent this problem; that is, where  , = max{0,  min − Δ(, )/} is the weight factor, Δ(, ) represents the distance between the centroids of element  and element ,  is the side length of the elements, and  min is the filter radius.

Algorithm.
The detailed topology optimization procedures for the design of a compliant mechanism using the above formulation and the modified additive hyperelasticity technique are described in Figure 2.With the proposed technique, the element instability in low-density elements can be completely overcome without sacrificing the accuracy of structural response analysis.Note that the main Update elemental constitutive matrix by equations ( 6) and ( 8 by equation ( 4) computational cost is the Newton-Raphson iterations for solving the nonlinear equilibrium equations.The MMA algorithm proposed by Svanberg [21] is used as an optimizer for updating design variables.The termination criterion of the optimization processes is the following: the maximal difference between the design variables of two adjacent iterations is less than 0.01.

Numerical Examples
The proposed topology optimization methodology is applied to the design of two compliant mechanisms undergoing large displacement.The solid material properties for the compliant mechanisms are given as follows: Young's modulus is  0 = 100 MPa and Poisson's ratio is ] = 0.3.The filter radius for design sensitivities is  min = 1.8.stiffness values  in = 1.0 N/mm and  out = 0.2 N/mm are put at point A and point B, respectively.The design problem is to find a force-inverting mechanism that maximizes the output displacement  out at point B. The material volume is restricted to 25% of the total volume of the design domain.

Design of an Inverter
In view of the symmetry, one half of the design domain is discretized by 5000 (100 × 50) 4-node plane elements.The   4(a)-4(c).The black region indicates the fully solid material, while the white region represents the low-density or void elements.It is noted that the optimal topologies change as the input force increases.For a small input  in = 2 N, the optimum material distribution is shown in Figure 4(a) with small output displacement of  out = 1.49mm.When  in = 36 N, the topology optimization procedure yields a rather different design as depicted in Figure 4(c).The main difference is the length and the incline angle of the right bars connected to the output point.In this case, the maximum output displacement is predicted to be  out = 23.15mm. Figure 5 plots the iteration histories of the proposed optimization processes.It can be seen that the optimization objective for each case converges to a stable value within 80 iterations.The corresponding details of solutions, such as the number of iterations, the maximum value of parameter   2 , the input displacement  in , the input work in ), the output displacement  out , the output work  out = (1/2) out  2 out , and the work transmission ratio  out / in , are listed in Table 1.It is observed that the topology in Figure 4(a) gains the most energy-efficient way to transfer input to output.However, it cannot generate large displacement at the output point.
The deformed shapes of the optimized topologies for three cases of input force are shown in Figures 6(a)-6(c).It is seen that some low-density elements in Figure 6(c) for the case of  in = 36 N undergo very large strains.However, these elements do not suffer from instability because the additive hyperelastic material is added.The final contour of the parameter   2 for the case of  in = 36N is plotted in Figure 7.In most parts of the design domain, the parameter   2 is extremely low.Therefore, the additive material only has slight effects on the overall response of the compliant mechanism.In order to show this, the optimal topology in Figure 4(c) was also analyzed by the nonlinear finite element method incorporating the element removal technique, in which the low-density elements with   ≤ 0.01 are removed from the finite element model, as shown in Figure 8.The curves of output displacement versus input work predicted by the modified additive material model are compared with the element removal model in Figure 9.In this largedisplacement case, the maximum discrepancy between the two sets of results is no more than 1.0%.On the contrary, if one adopted the previous method as suggested in [20] that uses a uniform value of  2 = 2.57 × 10 −6  0 for all elements, the maximum discrepancy becomes 4.7%.Therefore, the modified additive hyperelasticity technique is capable of eliminating the element instability and providing more accurate response analysis than the previous method during the topology optimization of large-displacement compliant mechanisms.

Design of a Gripper
Mechanism.The second benchmark example considered in this study is the topology optimization of a gripper mechanism.The geometrical dimension of the design domain and the boundary condition are given in Figure 10.The spring parameters are set as  in = 0.4 N/mm and  out = 0.1 N/mm.Here, the output displacement is to be maximized under two different input force values  in = 2 N and  in = 20 N. The allowable material volume ratio is specified as 25%.In the finite element model, only the symmetric upper design domain is discretized by 4839 nodes and 4688 four-node finite elements by using an element size of 1 mm.
The optimized solutions obtained by the proposed method for the two input force values are presented in   Figures 11(a) and 11(b).The resulting optimal solution in Figure 11(a) has a significant topological difference from the solution in Figure 11(b) due to the large-displacement effect.The latter design gains output displacement of 12.91 mm, while the former gains only 1.44 mm.The curves of output displacement versus input work for both designs are compared in Figure 12 based on the modified additive material model and the element removal model.The comparison shows that the relative discrepancy between the additive material-based response and the true response is rather small.In addition, even if a large input work ( in = 82 Nmm) is imposed, the design in Figure 11(a) can only yield 9.17 mm output displacement, which is much less than that of the design in Figure 11(b).This reveals that it is meaningful to account for the large-displacement effects in topology design of compliant mechanisms.
The deformed shape of the additive material-based model for the design in Figure 11(b) is shown in Figure 13, from which the effectiveness of using the proposed method to prevent instability in low-density elements can be validated.
It should be noted that, besides the element instability, another main difficulty in the topology optimization of compliant mechanisms is that there are some usually small hinge-like regions in the final topology, for example, in Figures 4 and 6.How to avoid such hinges during the topology optimization has been extensively investigated up to now.Many efficient methods, such as the length-scale control method [29], the morphology-based restriction scheme [30],  and the level set-based control method [31], have been proposed.Interested readers are referred to some recent publications [7,[32][33][34].

Conclusions
Based on the density penalization model and the modified additive hyperelasticity technique, a topology optimization method for large-displacement compliant mechanisms is presented.The modified additive hyperelastic material can effectively overcome the element instability difficulty, while exerting only a negligible influence on the structural nonlinear response prediction.Compared to the previous method, this modified technique provides more accurate results in structural response.The present method is easy to be implemented, without causing extra computational burden.It can generate a meaningful and well-defined topology when the mechanism is required to deliver relatively large displacement.Topology optimization of two benchmark examples of compliant mechanisms is tested.It is shown that the obtained optimal designs with large displacement have distinct differences in topology and performance compared with those with small displacement.

Figure 2 :
Figure 2: Flowchart of the topology optimization incorporating modified additive hyperelasticity technique.
Mechanism.The design domain and the boundary conditions of the displacement inverter are sketched in Figure 3.The top left corner and the bottom left corner of the design domain are fixed and an external force  in is applied in the middle (point A) of the left side.The domain size is 100 mm × 100 mm.Two artificial springs with

Figure 3 :
Figure 3: Design domain of the inverter mechanism.

Figure 5 :
Figure 5: Optimization histories for the inverter mechanism.

Figure 7 :
Figure 7: Contour of the parameter   2 for the inverter mechanism in the case of  in = 36 N.

Figure 8 :Figure 9 :
Figure 8: Element removal model for the design in Figure 4(c).

Figure 10 :
Figure 10: Design domain of the gripper mechanism.

Figure 11 :
Figure 11: Comparison of optimal layouts of the gripper mechanism with different inputs.(a)  in = 2 N. (b)  in = 20 N.

Table 1 :
Solution details for the inverter mechanism by the proposed method.