Robust design of large-displacement compliant mechanisms

. The aim of this article is to introduce a new topology optimisation formulation for optimal robust design of Micro Electro Mechanical Systems. Mesh independence in topology optimisation is most often ensured by using ﬁltering techniques, which result in transition grey regions di ﬃ cult to interpret in practical realisations. This problem has been alleviated recently by projection techniques, but these destroy the mesh independence introduced by the ﬁlters and result in single node connected hinges. Such features in the design are undesirable as they are not robust with respect to geometric manufacturing errors (such as under / over etching). They can be avoided by optimising for several design realisations which take into account the possible geometry errors. The design variations are modelled with the help of random variables. The proposed stochastic formulation for the design variations results in nearly black and white mechanism designs, robust with respect to uncertainties in the production process, i.e. without any hinges or small details which can create manufacturing di ﬃ culties.


Introduction
The focus in this article is on the design of compliant mechanisms by topology optimisation. Compliant mechanisms gain their mobility from the flexibility of the building components and they have found wide applicability in the production of Micro Electro Mechanical Systems (MEMS)small mechanical devices coupled with electronic circuits. The manufacturing is based on etching techniques utilised in the semi-conductor industry. The dimensions of MEMS are in the order of several hundred µm and due to their small size any hinges or assembly procedures are undesirable.
Topology optimisation (Bendsøe and Sigmund, 2004) has been utilised widely in the industry for optimising machine elements and assemblies. It is an iterative process where the aim is to minimise predefined objective, such as weight, cost or compliance, by distributing material in the design domain and fulfilling prescribed constraints. The design domain is discretised by using cells in 2D or voxels in 3D, and a design variable is assigned to each of them. The variables can take values 1 or 0, where 1 is assigned if a cell is filled Correspondence to: B. S. Lazarov (bsl@mek.dtu.dk) with material and 0 if it is void. In order to utilise gradient based optimisation methods, the 0/1 design problem is relaxed and the design variables are allowed to take values continuously between zero and one. The optimisation problem is mesh dependent 1 and convergence for mesh refinement is ensured by regularisation. Among the different methods proposed in the literature, the so-called filtering techniques have gained popularity. Initially, filtering has been introduced on the sensitivities of the objective (Sigmund, 1997), and later on the density field (Bruns and Tortorelli, 2001;Bourdin, 2001). The regularised topology optimisation problem results in designs with grey transition regions between the void and solid. These regions can often be removed by postprocessing, however in many cases they model the correct physics of the problem and discarding them will compromise the design performance. Recently, several projection schemes (Guest et al., 2004;Sigmund, 2007;Xu et al., 2010) have been proposed to decrease the grey transitions in the final designs. The first two (Guest et al., 2004;Sigmund, 2007) impose length scale on the void or the solid phase by regularising with finite support density filter and threshold projection with threshold 0 or 1, respectively. The projection scheme proposed in Xu et al. (2010) is based on Heaviside projection with arbitrary threshold. It results in nearly black and white designs with small features which are mesh dependent.
Topology optimised designs for MEMS consist of solid elements connected with hinges. The hinges in the case without Heaviside projection appear as grey material regions, and for projected designs as solid elements connected through a single node (e.g. Pedersen et al., 2001). The robust formulation (Sigmund, 2009;Wang et al., 2011), provide nearly black and white designs without any hinges by requiring the performance to be insensitive with respect to production errors in the geometry. The formulation is able to represent constant uniform under-or over-etching error distributed uniformly along the perimeter of the design. Under-or overetched design realisations are obtained using two different threshold projections. A more realistic representation of the manufacturing uncertainties requires to model them in a continuous space, i.e., considering the threshold to vary continuously between the most dilated and the most eroded case. The scenario can be modelled by using min/max formulation with more than three design realisations. Increasing the number of the design realisations will approximate closer the design space. A more systematic approach, presented here, is based on modelling geometry uncertainties by using stochastic theory. The threshold is modelled as a random variable. The linear elasticity state problem becomes non-deterministic and the methods developed for solving Stochastic Partial Differential Equations can be utilised to obtain the system response. The formulation proposed here is based on the stochastic moments of the mechanism response, and a full reconstruction of the system solution in the stochastic and the physical space is not necessary. The moments can be estimated easily by using Monte Carlo Simulations. The method converges relatively slowly to the true moments. As demonstrated in Xiu and Hesthaven (2005) an order of magnitude faster convergence for a limited number of the random dimensions can be obtained by the Stochastic Collocation Method, which is the solution method used for obtaining the presented results.
The paper is composed as follows. First, the standard deterministic topology optimisation approach is presented for large displacement linear elasticity in Sect. 2. The section covers the optimisation formulation, regularisation and projection techniques, and derivation of the objective sensitivities. In Sect. 3, the existing robust formulations are discussed and then the stochastic robust formulation is introduced. A brief discussion of the solution techniques for the stochastic state problem is presented in Sect. 4 and robust designs for MEMS obtained by the proposed approach are shown in Sect. 5.

Topology optimisation of large displacement compliant mechanisms
The objective in the considered compliant mechanism design, shown in Fig. 1, is to minimise the displacement in a selected degree of freedom. The original non-robust deterministic problem (Sigmund, 1997;Pedersen et al., 2001) can be written in discrete form as where the state elasticity problem is assumed to be discretised using the finite element method (FEM), u is the system response displacements vector, ρ is a vector with the topology optimisation variables associated with each finite element, and l is a vector with size equal to the size of u. The element l i which corresponds to the displacement degree of interest is set to one, and the rest are set to zero. r(ρ,u) is the residual vector function for the state problem. For a small displacement linear elasticity formulation of the state problem, the residual vector is given as r = f −Ku, where f is the external load and K is the stiffness matrix. For large displacement formulation r is presented in Sect. 2.1. The volume of the design domain occupied with material is denoted with V (ρ) and it is restricted to be smaller or equal to a predefined value V * . The design variables ρ i ,i ∈ N e are bounded between zero and one. The individual element contributions to the tangent matrix K are calculated by using elasticity modulus E obtained by the so-called solid isotropic material interpolation with penalisation (SIMP), which can be written as where E 0 is the stiffness of the solid phase. E min is the stiffness of the void phase -a small number larger than zero in order to ensure non-singularity of the tangent matrix. The parameter p is used for penalising intermediate design values and is usually taken to be p = 3, andρ i in Eq. (3) is the physical density at the selected point in the design domain. In order to ensure mesh independence of the optimised solution, as well as to avoid checker-boards, the original design field ρ is filtered. The filtered densityρ can be obtained explicitly by using weighted average of the design variables around each element (Bruns and Tortorelli, 2001;Bourdin, 2001) or implicitly by solving partial differential equation (PDE) for the filtered density field  − r 2 ∇ 2ρ +ρ = ρ, ∂ρ ∂n = 0 where r is a filter length parameter which defines the length scale imposed by the filter. For r = 0 the filtered field is equal to the original design field. The increase of the filter parameter r suppresses fast oscillations in the design field and passes out the slowly varying components of the field. The boundary condition ensures that the filter is volume preserving, i.e. the volume of the input design field ρ is equal to the volume of the filtered fieldρ. The vector n denotes the outward normal to the design boundary. The PDE (4) is discretised by using the same mesh utilised for solving the state problem, and in discrete form can be written as where K f is the discrete differential operator ∇ 2 + 1 and T f maps the design field vector associated with each element to the nodal input field of the PDE filter. An example of a MATLAB implementation of the PDE filter can be found in Andreassen et al. (2011). If the physical element density is represented by the filtered density obtained by solving Eq. (4), the optimised design consists of grey areas which are difficult to interpret. Practical realisations require a discrete black and white solution. Such a solution can be obtained by threshold projection (Guest et al., 2004;Sigmund, 2007;Xu et al., 2010). All values above the selected threshold η are projected to 1 and all values below the threshold are projected to 0. Mathematically the operation can be represented by the Heaviside function, which is not differentiable. Therefore, for computational purposes, it is replaced by a smooth function with the expression suggested in Wang et al. (2011) and given aŝ In the limit when, β → ∞ Eq. (6) approaches the Heaviside function with threshold η. One undesirable feature of the threshold projection is that the length scale imposed by the density filter is lost. In Wang et al. (2011), this property is demonstrated for several optimisation problems in heat transfer and compliant mechanisms designs. The optimised designs consist of small features comparable with the mesh size. Furthermore, the threshold projection for the compliant mechanisms results in hinges in the final black and white design which is not desirable.

Non-linear elasticity and finite element formulation
It is assumed that the mechanism displacements are large and the standard small displacements and small strains formulation in linear elasticity is not capable of representing the final deformed state of the system. In order to account for finite deformations of a continuous body, the linear strain and the Cauchy stresses are replaced with a non-linear strain measure and its conjugate stress. Detailed overview and finite element discretization for finite strains elasticity can be found in many textbooks on the subject (e.g. Krenk, 2009;Belytschko et al., 2000;Bonet and Wood, 1997). Here the non-linear strains are considered to be the Green's strains given in a tensor form as where D is the displacement gradient tensor with respect to the initial coordinate system. Each component of E can be written as where Einstein summation convention is assumed with respect to the index γ. By removing the second quadratic term in Eq. (7), the small strain measure is recovered. Green's strains require the introduction of the 2nd Piola Kirchhoff stress tensor. Both of them are work conjugate. A linear constitutive relation is assumed between the strain dE αβ and the stress dS αβ increments in the form where the material tensor C is obtained as with E given by Eq. (3), and C 0 αβγδ -the material tensor for unit elasticity modulus.
An expression for the residual forces is obtained by forming the virtual work equation and requiring that the variation of the total work is zero p is an external force vector obtained by integrating any point, surface or volume forces acting on the system. The stress vector s consists of the following stress tensor components s = [S 11 ,S 22 ,S 33 ,S 23 ,S 13 ,S 12 ] T (12) B(u) is a matrix function which depends on the current deformed state and relates the strain vector variations to the displacement variations In total Lagrangian formulation the integration in Eq. (11) is performed over the original undeformed volume. At equilibrium the residual vector is equal to zero, and the solution of with respect to u determines the deformed state of the system. The values of u are obtained iteratively by the Newton-Raphson method with tangent matrix computed as www.mech-sci.net/2/175/2011/ Mech. Sci., 2, 175-182, 2011

Optimisation sensitivities
The objective sensitivities can be obtained by using adjoint analysis, and for large displacement formulation detailed derivations can be found in Pedersen et al. (2001). The gradient of the objective with respect to the physical design field is given as where λ is obtained as a solution of the following system of linear equations The matrix K t corresponds to the tangent matrix computed at the state equilibrium, i.e. for r(u) = 0, and l is an input to the system which is zero everywhere except at the degree of freedom where the objective is computed. The above derivation is based on the assumption of symmetry in the tangent matrix.
The gradients with respect to the design variables ρ are computed using the chain rule The derivative ∂r/∂ρ i is computed analytically by differentiating Eq. (11) with respect toρ. With respect to the filtered field it is computed by applying the chain rule and after that analytically differentiating the threshold projection given by Eq. (6). In discrete vector form the gradients with respect to the nodal values of the filtered field are given as where s is assembled element-wise by integrating the sensitivities contribution from each element . The final gradients with respect to the original design variables associated with the elements are computed as

Robust topology optimisation
A general overview of various formulations for obtaining robust solution to an optimisation problem can be found in Beyer and Sendhoff (2007); Tsompanakis et al. (2008). In this work, robustness is required for the system performance with respect to uniform under-or over-etching of the design, i.e. the optimised design has to perform well when the mechanism elements are produced thinner or thicker with respect to a reference topology supplied to the manufacturer. The geometric variations in the design topology can be modelled by varying the threshold η in Eq. (6). Three projections with three different thresholds η e ,η i ,η d are shown in Fig. 4. The three projections are called eroded, intermediate and dilated (Sigmund, 2009;Wang et al., 2011). If the intermediate projection is considered to be the reference design, uniform over-etching error can be modelled by the difference between the intermediate and the eroded design projections. Uniform under-etching error can be modelled by the difference between the intermediate and the dilated design projections. Using these three cases, robust designs for small displacement compliant mechanisms are obtained in Wang et al. (2011) by minimising the maximal objective of the three projections. The formulation is an extension of an earlier work by Sigmund (2009) where the min/max formulation is applied for eroded and dilated designs obtained with thresholds η = 1 and η = 0, respectively. Both formulations utilise several discrete points in the design space and they do not account for the continuous nature of the geometric errors. A way to expand the considered error space is to use the min/max formulation for more than 3 cases, however the latter would increase the computational burden significantly.
A continuous geometric error can be modelled systematically by employing a stochastic variable with suitable (physically admissible) distribution function. Here it is modelled by using the threshold projection Eq. (6) where the threshold η is considered to be uniformly distributed, i.e. η ∈ U[a,b]. The lower bound a of the uniform distribution corresponds to the most dilated case, and the upper bound b -to the most eroded case. The mean threshold corresponds to the reference design supplied to the manufacturer. Representing the threshold as a random variable results in random variations in the deterministic objective considered in Eq. (1), as well as in the volume occupied with material. The original deterministic optimisation Eq. (1) can be reformulated by using a probability measure of the mechanism performance, or the moments of the objective distribution. By utilising the stochastic moments of the response, the following stochastic robust optimisation problem can be introduced where E l T u is the expected value of the deterministic objective and STD l T u is the standard deviation. Alternatively optimal robust design can be obtained by The first optimisation problem given by Eq. (22) minimises the mean performance by constraining its standard deviation.
The second formulation Eq. (23) provides an alternative where the parameter κ controls the contribution of the standard deviation to the objective. Increasing κ puts more weight on the standard deviation and the obtained solution possesses response which is less sensitive to geometry variations. In both formulations, based on the results from Sigmund (2009) and Wang et al. (2011), the volume constraint is imposed on the dilated design, i.e. V d (ρ) ≤ V * . The dilated design has the largest amount of material in the considered model.
The expectation E l T u of the original deterministic objective can be computed as where ϕ(η) = 1/(b − a) is the probability density function for uniform distribution U[a,b]. The standard deviation of the deterministic objective is given as STD l T u = Var l T u , and the variance is computed as

Optimisation algorithm and numerical implementation
The main difficulty from computational point of view, in the stochastic robust formulation, is the evaluation of the mean and the variance of the deterministic objective. These can be estimated by obtaining a solution of the finite strain elasticity problem with stochastic modulus of elasticity. Several solution strategies can be employed (Xiu, 2010), and among them the easiest and the most expensive one for a single random variable, in terms of computations, is the Monte Carlo simulations (MCS) method. The method converges to the true expected value relatively slow, with a rate proportional to the inverse of the square root 1/ √ M of the number of the realisations M. For a sufficiently smooth solution of the stochastic partial differential equation problem, the Stochastic Collocation Method (SCM) (Xiu and Hesthaven, 2005;Xiu, 2010) converges an order of magnitude faster than MCS. The method doesn't require the actual construction of the solution in the stochastic space, and the integrals for the expectation Eq. (24) and the variance Eq. (25) can be computed by evaluating the solution of the deterministic problems at prescribed collocation points η i . The SCM method is based on Lagrangian polynomial approximation in the stochastic space. The residual of the interpolated solution is required to be zero at selected collocation points, which can be written as Each one of the above equations is equivalent to the state problem formulated for threshold η k . The integrals Eq. (24) Methods in Applied Mechanics and Engineering, 190, 3443-3459, 2001. Guest, J., Prevost, J., and Belytschko, T.: Achieving minimum length scale in topology optimization using nodal design variables and projection functions, International Journal for Numerical Methods in Engineering, 61, 238-254, 2004. Krenk, S.: Non-linear   and Eq. (25) can be computed as where ω k are integration weights. The sensitivities for the stochastic robust formulation Eq. (22) where c k = c(ρ,η k ). The gradients ∂c k /∂ρ i can be estimated using the derivation presented in Sect. 2.2.

Numerical examples
The proposed stochastic robust formulation is demonstrated for the design of compliant inverter mechanism, with design domain and boundary conditions shown in Fig. 1. The length of the design domain is chosen to be L = 300 µm and the thickness is t = 7 µm. The elasticity modulus is set to E 0 = 180 GPa. The input and the output springs stiffness is set to k in = 4.00 mN µm −1 and k out = 0.01 mN µm −1 , respectively. The driving force is f in = 20 mN. The integration of the expectation and the variance is performed using Gaussian quadrature. The error tolerance for the Newton-Raphson iterations is set to be 10 −6 and the algorithm is stabilised using arc-length control. The Method of Moving Asymptotes (MMA) (Svanberg, 1987) is utilised for solving the optimisation problem. All designs are obtained using a continuation scheme with respect to the projection parameter β (see Eq. 6). The optimisations start with β = 1 and every 50 steps, β is increased with 1. When β is equal to 16, the continuation scheme doubles it and the design process runs for 100 iterations. In order to decrease the computational cost, all steps except the last one (β = 32) are performed with 2 Gaussian integration points. Such a numerical integration scheme is not able to approximate the standard deviation well, however based on numerical simulations, it captures very well the mean response. The initial steps are used only to obtain good initial guess for the final optimisation step, thus it is not necessary to estimate the moments of the response very precisely. The final step of the optimisation (β = 32) is performed with Kronod-Patterson quadrature using 31 points. The optimisation process can be improved further by using the nested property of the Kronod-Patterson integration points and estimating the error in the integration process. The robust designs are relatively stiff compared to the non-robust and small displacement theory can be used in the initial steps to decrease further the computational cost. The difference in the performance for designs obtained by using large and small displacement formulation, in the examples presented here is between 15 % and 20 %. Therefore for small β linear analysis can save significant amount of computational time. For large β the design changes relatively slow and if the initial guess if far from the optimal one, the optimisation would require large number of iteration steps. In order to avoid half width elements close to the design domain borders (e.g., Wang et al., 2011), the boundary conditions (BC) for the PDE filter (Fig. 2) differ from the original formulation . Dirichlet BCρ = 0 is imposed on two sides of the design domain andρ = 1 is imposed around the input and the output of the system. The BCρ = 1 implies solid material outside of the design domain andρ = 0 implies void.
The first example is an optimised topology of the compliant mechanism using the large displacement deterministic formulation. The result is shown in the middle of Fig. 3, and is obtained with threshold projection η = 0.5. In order to use the same settings as the ones for the robust formulation, the volume constraint is imposed on design obtained with threshold projection η = 0.3. The optimised topology consists of elements and one node connected hinges. Small deviations from the design are shown on the figure as well. The left design in Fig. 3 is obtained by erosion with Heaviside projection threshold η = 0.6. The mechanism is completely disintegrated. The right design in Fig. 3, is obtained by dilation with projection threshold η = 0.4. The mechanism hinges are filled with material and the mechanism becomes very stiff compared to the reference one obtained for η = 0.5. The mean value of the deterministic objective for η ∈∈ U Three projections for an optimised design obtained by using large displacement robust formulation are shown in Fig. 4 and Fig. 5. All of them perform similar for the selected thresholds. The mean for η ∈ [0.4,0.6] is E[c] = −9.58 and the standard deviation is STD[c] = 3.46. Furthermore, in contrast to the non-robust deterministic design, the mechanism does not posses any hinges, and small erosion or dilation does not disintegrate it. The penalty is smaller maximal displacement. The performance and the standard deviation of the design for η ∈ U[0.3,0.7] is E[c] = −8.96 and STD[c] = 5.80. The mean performance in the design interval is better for the robust design compared to the one for the deterministic case. In addition it is robust with respect to erosion or dilation, i.e. the standard deviation is smaller for the design obtained by using the robust formulation. The projections shown in Fig. 4 are obtained with two integration points and the projections in Fig. 5 are obtained with 31. The topology in the second case (Fig. 5) differs slightly from the one obtained with two integration points. The performance is improved slightly and the main difference is in the eroded design which is thicker than the one shown in Fig. 4. This behaviour is expected, as more precise approximation is based on integration points which are closer to the bounds of the threshold interval and therefore the sensitivities for highly eroded or highly dilated structures will have contributions to the average sensitivities given by Eq. (29) and Eq. (30). Threshold projections for a design obtained with broader threshold interval η ∈ [0.2,0.8] are shown in Fig. 6. Increasing the threshold interval decreases the performance of the mechanism and increases its robustness with respect to geometry variations. The optimal threshold interval, as well as the selected threshold distribution, have to be tuned to a given production process. Increasing V * makes the volume constraint inactive. For the selected formulation and boundary conditions, the optimisation finds topology with the highest possible flexibility rather than the highest possible force transfer.                     www.mech-sci.net/2/175/2011/ Mech. Sci., 2, 175-182, 2011

Conclusions
A new optimisation procedure for robust design of compliant mechanisms is demonstrated. The performance of the obtained designs is robust with respect to uncertainties in the geometry. The uncertainties are modelled using Heaviside projection with a random threshold which is selected to be uniformly distributed in the threshold interval. The obtained designs do not possess any hinges and the requirement for robustness ensure easy manifacturability. For large complex models the proposed formulation needs further improvements in order to decrease the number of optimisation iterations and the computational cost associated with each of them. This will be subject of future work.