Stress-based topology optimization of compliant mechanisms using nonlinear mechanics

. The present work demonstrates how a light structure can be easily designed through Topology Optimization even including complex analysis and sizing criteria such as hyperelastic Neo-Hookean materials for nonlinear analysis and aggregated stress constraints. The SIMP approach was adopted and two different strategies were analysed using an in house versatile MATLAB code. MMA was used as reference optimizer (in structural optimization) whereas a unified aggregation and relaxation method was adopted to deal with stress constraints. Feasibility was analyzed from the viewpoint of allowable stress verification. Two test cases are then studied: a morphing airfoil (for aeronautical applications) and a geometric inverter (for mechanics and bio-medical applications). For both, a hyperelastic Neo-Hookean material was chosen. Finally a complementary study on the effects of constraints and the input force intensity is also presented.


Introduction
The past twenty years have seen increasingly rapid advances in the field of structural optimization: in particular, Topology Optimization (TO) is becoming a key figure in the research panorama.This method leads to an optimal structural efficiency, through the removal of unnecessarily placed material, in order to accomplish the minimization of a given objective (such as the compliance of the system).Since the pioneering work of Bendsoe and Kikuchi (1988) [1] a number of approaches have been proposed, among which are: density-based methods (the well-known SIMP algorithm) [2,3], evolutionary strategies [4] and level-set based methods [5,6].
Moreover, several researchers have tried to implement and diffuse simple and efficient codes.One of the most well-known and important results is due to Andreassen et al. [7], who created a very simple and efficient 88-line code written in MATLAB, capable of optimizing an MBB beam through a SIMP algorithm.Compared to other available codes and the same 99-line code (published by Sigmund et al. [8]), the former needs much less computational power, thanks to judicious use of vectorization and pre-allocation of most of the variables.The present work * e-mail: gabriele.capasso@student.isae-supaero.frfocuses on the SIMP approach and the code developed in our study is largely inspired by the top88.
In the literature, several promising applications of TO to classical aerospace structures can be found [9][10][11][12][13].Recently, researchers have shown an increased interest in alternative architectures presenting wings with a variable geometry, better known as "morphing wings" [14][15][16][17].Different types of morphing have been tested both numerically [18] and experimentally [19], revealing the advantages in terms of aerodynamic drag, noise and vibration control with respect to conventional aircraft configuration, based on flight control surfaces.Two main categories can be distinguished [20]: planform morphing, based on the modification of the plant itself of the wing, and performance morphing, which leads to variating the airfoil in shape (shape morphing) or in incidence (twist morphing).In this paper, a shape morphing design is analyzed.
Earlier work on TO was restricted to linear structural designs.Recent research explored the influence of nonlinear mechanics on the TO in several simple applications, including geometrical nonlinearity [21][22][23][24], material nonlinearity [25][26][27][28][29][30] and both geometrical and material nonlinearities simultaneously [31][32][33].Work from [24] evidences the differences encountered in a TO problem, according to the choice of a linear or a nonlinear analysis, leading to totally different results.The introduction of Nonlinear Mechanics into the TO of a morphing wing has already been treated by Bhattacharyya et al. [9]: in particular, geometric and material nonlinearities have been investigated.In [9], the optimization problem presents the objective of minimizing the actuating force, under constraints of desired volume fraction, compliance and trailing edge displacement.The obtained structure is, in fact, a compliant mechanism capable of getting a geometric advantage.The first of the two parts of the present work focuses on the same structure, whilst the objective is to design the lightest structure under stress constraints, in order to conceive a feasible structure.
In recent years other compliant mechanisms have frequently been proposed for aerospace systems due to their advantages over traditional rigid-link mechanisms.As compliant mechanisms can be manufactured using fewer parts (usually one single part), they are lighter than their traditional counterparts.They are also easy to manufacture with modern 3D printing techniques, and do not require assembly.The ability to 3D print mechanisms and tools on board the International Space Station (or future bases on the Moon and Mars) will save substantially on the cost and time involved in transporting these devices from Earth.To be prepared for deep-space exploration and habitation in the future, research and technology development concerning the utilization of additively construct structures from insitu materials is incorporated into the in-space manufacturing technology roadmap [34].
A compliant mechanism has been investigated by Conlan-Smith and James [35].Their study focuses on the optimal design via Topology Optimization of a bioinspired structure.The second part of the present work is inspired from [35], and aims at solving the same optimization problem by considering feasibility constraints, in terms of admissible stresses.
Stress-based Topology Optimization is an active research field.Several examples available in literature show the challenge posed by stress constraint in TO problems.Essentially, two problems arise [36]: first, a large number of constraints must be considered, since unlike stiffness, stress is a local quantity; second, stress is highly nonlinear with respect to design variables.Recent advances show possibilities to overcome these difficulties [37][38][39][40][41][42].There are few recent examples in the literature where the stress distribution is considered as a constraint in TO problems with geometrical nonlinearities (as for instance [39]).The introduction of both stress constraint and nonlinear mechanics into a TO using nonlinear mechanics is the main challenge of the present work.
The present work aims at designing two compliant mechanisms: a morphing wing and a geometric inverter.The former is treated as an academic example conceived to validate the method, while the latter is a more realistic test-case.In both mentioned parts of this article, the optimal design is achieved using a nonlinear analysis of a hyperelastic material through TO, considering the stress as a constraint of the problem.Section 2 illustrates the theoretical background behind the optimization problem.Section 3 details the mathematical formulation and the newly introduced elements of the paper.Section 4 presents the different test-cases, outlining the relative mathematical formalization.Optimized designs and numerical results are finally reported and discussed in Section 5.

Theoretical background
In this section the basic principles of nonlinear mechanics and Topology Optimization are explored.

Nonlinear mechanics
The present work focuses on the design of largedisplacement mechanisms, in which material doesn't follow traditional linear mechanics.This topic will be just introduced in this section, although it is largely explained in [43].
The resolution of a Nonlinear Finite Element Analysis (NFEA) is not straight-forward, since the global system is nonlinear: in [43] several methods are illustrated.The most common algorithm is the one which goes under the name of Newton-Raphson, which is usually modified through the introduction of a step-force control.This method guarantees better stability than the simple Newton-Raphson, preserving the quadratic convergence [43].
The former is the main subject of the present article and will be explained in the following section.

SIMP algorithm
The density-based approach, also known as Solid Isotropic Material with Penalization (SIMP), is the most popular TO method.SIMP has a solid mathematical foundation [1,[47][48][49][50].This method is capable of handling various objectives and constraints, and is relatively easy to implement within a finite element environment.
The basic idea behind NFEA based SIMP is that each finite element is associated with a fictitious pseudo-density variable 0 ≤ ρ ≤ 1, that essentially parameterizes the topology.The pseudo-densities are optimized to reach the desired objective.The general algorithm is detailed in Algorithm 1.In [51] the reader can find an exhaustive review on numerical instabilities typical of TO.They include essentially checkerboard patterns, mesh-dependencies and local minima.To deal with the first two issues, a meshindependency filter can be applied [51].A fundamental aspect of such filters is the distinction between the physical densities ρ and the merely numerical ones x: the latter category constitutes the design variable array updated in the optimization process at each iteration; once having updated numerical densities, these are filtered to obtain the physical distribution of the material.NFEA is based on physical densities ρ.

Optimizer: method of moving asympthotes
When the number of design variables is very large, gradient based optimization methods are the most efficient algorithms.They are preferred to gradient-free optimization methods since the latter typically require many evaluations.
As shown in [9], the MMA is one of the most suitable optimizers which could substitute the Optimality Criteria adopted in the 88-line code by Sigmund [52], in order to accelerate convergence.
This method is based on a local convex approximation of objective and constraint function: this needs as input the local evaluation of the objective function, the constraints and all their derivatives with respect to the numerical densities x.
The solution is said to have reached convergence when the first order Karush Kuhn Tucker (KKT) conditions are satisfied to an absolute tolerance of 2 × 10 −3 .

Material description
A nonlinear hyperelastic material was investigated: the same methods described in [9] were adopted, but modifying the proposed algorithm in order to achieve the same results as well as a fast convergence.Moreover, far too little attention has been paid to the feasibility of the structure itself: this can be observed by the absence of a constraint related to the maximum stress supportable by the material.
The chosen material is described by the potential density defined in [9,53] and reported in the formula below: Here J is the determinant of the deformation gradient F ; the term tr(C) denotes the trace of the right Cauchy-Green deformation tensor C = F T F ; finally, constitute the Lamé constants of the material which are seen as functions of Poisson ratio ν and Young Modulus E. This latter is supposed to be related to the density ρ of each element (e), through one of the two following relations: which will be used in order to make comparisons in our forthcoming NFEA.Following [7], equation ( 2) is referred to as calssical SIMP and equation ( 3) as modified SIMP.

Nonlinear finite elements analysis
We treat the NFEA as another minimization sub-problem in the global optimization: the function to minimize is the total energy of the system, given by both internal and external forces, while the constraints are the displacements imposed in the boundary conditions.The domain is decomposed using triangular elements.The definition of the objective function derives from the Principle of Virtual Works, where the external work is computed as the product between the force and the displacement of the point of application: Here, U represents the vector of the nodal displacements which has to be determined by Variational Principle.In order to make the Optimization Toolbox in MATLAB more efficient, we have to provide it with the gradient and the Hessian of the total energy.The gradient will be the residual R, computed as the difference between internal and external forces: where The Hessian will be simply the tangent stiffness matrix K T , computed as the second derivative of the compliance to the nodal displacements: The equilibrium condition is satisfied when: The MATLAB tool fmincon is a suitable instrument to perform NFEA, since it treats this sub-problem as a constrained minimization problem.Moreover, it contains a stable implementation of the Newton-Raphson Method in the Interior-Point Region [54].This feature turns the fmincon process into a fast and reliable solver for nonlinear mechanical problems.

Stress constraint treatment
The Second Piola-Kirchoff stress S relates forces in the reference configuration with areas also measured in the reference configuration (see [53]).It can be calculated as twice the derivative of the potential energy function with respect to the right Cauchy-Green deformation tensor: The Cauchy stress, σ, which relates the forces in the deformed configuration with areas in the deformed configuration, can then be defined as: As a consequence, the microscopic (local) stress tensor will be evaluated as: where B is the left Cauchy-Green tensor, which is defined as B = F F T , and I is the second-order identity tensor.
Here, microscopic stress was considered by calculating the Von-Mises stress, under the hypothesis of plain strain.It is important to underline the fact that this does not depend directly on density distribution.In this work the unified approach proposed by Verbart et al. [40] was adopted to incorporate stress constraints in topology optimization.
Assuming that the densities in SIMP represent a porous micro-structure, one can distinguish the stress at macroscopic (effective) versus microscopic levels [55].The former is the stress computed using the Young modulus of the SIMP model (c.f.Eqs. ( 2)-( 3).The latter is computed considering a stress model that mimics the behavior of the "local stress" in a rank-2 layered composite [55].For this reason we identify the microscopic stress as: where S 0 is the second Piola-Kirchoff stress considering a material with full Young Modulus E = E 0 (i.e.ρ = 1).

Adjoint sensitivity analysis
The method which is going to be illustrated has been proposed by Bhattacharyya [9].Since the objective and constraint function depend both on density distribution and nodal displacements, an adjoint sensitivity analysis is to be performed.In particular, we use an augmented Lagrangian L[ρ, U F (ρ)], which may be defined as follows: Here f indicates a general function, indifferently the objective or the constraint; λ is a Lagrange multiplier ensuring that the residual of forces R vanishes as required in equation (8).Deriving this expression and using the chain rule, we obtain: Collecting all the implicit terms indicated by d dρ , we obtain: In order to eliminate the implicit dependence of free nodal displacements on density distribution, the lagrangian multiplier can be computed, remembering the definition of the derivative of the residual related to the nodal displacements, as follows: Finally, we can obtain: One can observe that the evaluation of gradients only requires a linear system of equation resolution, extremely convenient if compared with direct or finite difference approaches.

Density filtering
As pointed out in Section 2.2, a mesh-independency filter was applied to deal with certain numerical instabilities, namely checkerboard patterns and mesh-dependent solutions [51].
Here the method suggested by Bruns and Tortorelli [56] was implemented, where: Here w ij is the weight associated with the ith element within the prescribed neighborhood of the jth one, determined as: The parameter R min , called filter radius, is the radius of the specified neighborhood and r ij the distance between the centroids of the two elements.This formulation proposed in [56] is valid under the assumption of uniform mesh.This filter makes the results independent from the adopted mesh, improving the resolution as R min decreases.However, the maximal resolution is linked to the mesh itself, since the minimum applicable value R min cannot be smaller than the size of the largest element.
When computing derivatives with respect to the numerical densities x (necessary to gradient-based optimization methods), two steps are followed: firstly, derivatives with respect to physical densities ρ have to be computed; then, chain rule is applied to obtain the desired gradient.

Constraint aggregation
Treating the relaxed stress distribution, there would be an excess of constraints.So the maximum may be taken into account to represent the whole domain.Given the fact that the max function is non-derivable, the lower bound Kreisselmeier-Steinhauser function [36,57] was employed: The larger the factor P , the better the approximation of the max, but also the more expensive the computational cost of the overall optimization process.In fact the non-linearity of the optimization problem is enhanced by the use of higher values of P .Therefore, P has to be chosen as an appropriate compromise between computational burden and accuracy in stress control.

First optimization problem: NACA airfoil
As an academic example, this part focuses on the minimization of the weight of a shape morphing wing.Using the classical SIMP approach (see Eq. ( 2)), the whole problem can be formalized as follows: Here V represents the volume fraction, U e the displacement at the trailing edge and σ lim is the highest tolerable stress in the structure (to avoid failure).The constraint relative to the σ in each element is due to the desire to design a feasible structure while not reaching the weakest one: in particular, we will consider the Von-Mises stress, under the hypothesis of plain strain.Finally, the constraint on U e derives from the necessity of designing an efficient shape morphing wing.We shall also perform further analyses based on the modified SIMP approach (see Eq. ( 3)), for comparisons.The new optimization problem is then: Here another constraint has to be added, since the volume fraction minimization under stress constraint, would generate gray regions, preventing the creation of a continuous chain between force application point and trailing edge.Such constraint will be an inferior limit on the compliance C of the trailing edge region [35]: this value is arbitrarily defined as the initial compliance.

Second optimization problem: geometric inverter
The second optimization problem is a more realistic test case, inspired by the work in [35].Figure 2 illustrates the initial domain.
The idea is to design a compliant mechanism which would invert the input displacement, creating as output another displacement on the opposite side of the structure.The actual input is a force.The optimized mechanism is defined as the one capable of reaching the highest ratio between output and input; as the two displacements have opposite signs, this is translated into a minimization problem: The second pseudo-density parameterization of the SIMP Algorithm is adopted.Before the main loop, the topology is initialized by imposing ρ = V 0 in every element of the domain.The constraints are quite similar to the ones adopted in the previous section.We have thus: -volume fraction: the structure has to be as light as possible, to reduce costs and weight; -stress: the whole mechanism has to support the applied load without breaking, or reaching the yield point; -compliance: a continuous chain between input and output is ensured, by imposing that the output region has to keep the same compliance at the first iteration loop (where ρ = V 0 in every finite element) [35].
Another aspect which appeared fundamental is the symmetric nature of the domain: to accelerate computation, only half of the whole domain has been simulated.

Results and discussion
In this section, numerical results are presented and discussed.Simulations have been performed using a unique MATLAB code, adopting the original MMA code provided by Svanberg [52] and without altering its default parameters.Given the fact that the constraint is approximated (as explained in Sects.3.3 and 3.6) stress distribution post-processing on the optimal structure is necessary.The advantage of our methodology is that stress distribution is already computed in our optimization loop ("Constraints evaluation" in Algorithm 1).

NACA airfoil
A first analysis is performed on the two problems 22 and 23, considering the same domain (represented in Fig. 1).The parameters adopted in this section are reported in Table 1.

First formulation results
In this subsection, the results to the problem (22) are presented.The base structure can be observed in Figure 3a.This can be decomposed into four main sectors: -Front sector: This is responsible for stress redistribution, as illustrated in Figure 3d.The stresses use this path in order to create an equilibrating moment.The middle part, situated between the front chain and the application point of the external force, is empty: in fact, the volume fraction constrained minimization produces this result, where the stresses have to be collocated in an intermediate location.Moreover, this sector connects the two series of points where the homogeneous Boundary Conditions are applied.-Central chain: This sector is the transmitter of forces between the application point of the external force and the rear chain (discussed later), as shown in the macroscopic stress distribution (Fig. 3d).It is directly linked to all the other parts of the optimized structure.-Inferior chain: It is the main link between the central chain and the inferior BCs.In the base structure, It is also directly linked to the front part: this is due to the fact that it starts exactly at the left limit of the inferior boundary points affected by the application of homogeneous BCs.-Rear chain: This is the only part not interested in any redistribution of stresses, as shown in Figure 3c  and d.Its main role is transmitting the commanded displacement at the trailing edge, thanks to an input force provided by the central chain.This sector affects the volume fraction negatively, given that it doesn't have to carry any stress.
The optimized structure in deformed configuration is illustrated in Figure 3b.The reader can observe the logic at the base of the mechanism.By introducing a single input force at the node in green, the effort is transmitted to the rear chain through the central sector.Simultaneously, the inferior chain and the front sector connect the central chain to the nodes in blue (where homogeneous Dirichlet Boundary Conditions are applied).Finally, the rear chain converts the input force into a displacement at the trailing edge.
The distribution of the relaxed constraint proves the approximated nature of the G KS function: this is not always guaranteed, but it is the product of a compromise between precision, stability and computational cost.However, the only elements where the relaxed constraints are not satisfied are, effectively, "empty".

Second formulation results
When using the second formulation to parameterize the topology, a totally different structure is obtained.Results from problem (23) are presented in this subsection.
The differences are evident by comparing Figure 3 and Figure 4.
As was visible in the previous numerical experiment, the aggregation method does not always guarantee the respect of the stress limits: this is a consequence of the G KS approximation.This time, an overshoot can be found on the filled structure.A lower limit on the stress should have been introduced in order to adapt to the real constraint.
The logic behind the optimal structure presented in Figure 4 (solution of the problem ( 22)) is unaltered compared with to the first one in Figure 3 (solution of the problem ( 23)).
However, this formulation leads to a more efficient application of the stress constraint, since all empty elements do not contribute to the macroscopic stress distribution.An immediate consideration is the removal of most of the Front sector and the Rear chain from the previous results: globally the obtained volume fraction would be inferior.
As a consequence, the modified SIMP approach (see Eq. ( 3)) is more suitable for topology optimization problems presenting a stress constraint.This consideration is applied in the second part of the present work.

Geometric inverter
Out of the previous academic examples, the second formulation has been identified as a better compromise between severity of the penalization on stress constraint and reality of the final structure, leading to a much lighter final structure.However, a lower limit on the stress has to be introduced.In the second test case of the present article, modified SIMP was adopted.Moreover, only half of the domain has been considered in the simulation to accelerate computation (symmetric properties have been exploited).All the parameters characterizing the Inverter problem are reported in Table 2.
In Figure 5 the final structure is presented.In the deformed configuration (Fig. 5b), the gain in terms of displacement inversion is easily visible.
Each half of the compliant mechanism can be decomposed into four main sectors: The relaxed stress constraint and macroscopic stress distributions present some differences: in fact, the former reports a stress constraint violation in the elements of the terminal chain of the actuator, where the displacements are higher, due to the aggregation; this violation disappears in the macroscopic stress distribution.On the contrary, the maximum effective stress is in the areas of higher curvature.

Influence of stress constraint
Further numerical analyses were performed, in order to investigate the effects of the stress constraint.In particular, both increases and decreases on the stress limits were considered.
If the stress limit decreases to 0.2 MPa, the obtained structure is totally altered (results reported in Fig. 6).However, the global decomposition follows the same logic adopted in the previous section.Four main differences may be noticed: -the Central Bulb is more compact and does not present any hole; -the Support Net presents only two ramifications; -the Angular Chain and the Support Net cannot be distinguished; -the Terminal Chain is much longer than before.
The reader can easily understand the necessity of introducing an additional buckling constraint if the stress limit is exaggeratedly low.
An analogue analysis was performed for an augmented allowed limit stress of the material (results reported in Fig. 7).The differences are less evident than an analogue reduction on stress limits.However, the Support Net is more ramified, allowing for a different allocation of the efforts.
Another fundamental aspect is the imperfect exploitation of the material: in fact, the stress constraint isn't active, leading to a structure where the maximum macroscopic stress isn't reached (refer to Fig. 7d).In Table 3 the numerical results are shown for each value of σ lim .In particular, it has been demonstrated that higher limit in terms of allowed stress doesn't necessarily imply an augmentation in final performances.Moreover, the number of iterations to convergence augments as σ lim decreases.

Influence of the amplitude of the applied load
The applied load shows little differences with respect to the base structure.In particular, the Support Net is less ramified for higher loads : this reflects again the same effect of a very small reduction in acceptable limit stress.
Moreover, as the amplitude of the force increases, an evident overshoot on the macroscopic stress can be  observed: this leads to the idea that another approximation for the max is needed.
In the Appendix, different designs of the optimized structure are presented.
In Table 4 the numerical results are presented for each value of F app .By increasing the amplitude of the applied load, the geometric advantage improves: however, the gain isn't proportional to the relative augmentation in terms of applied force.On the other hand, a general rule governing the number of iterations cannot be exploited.

Conclusion
Topology Optimization was used to design two different compliant mechanisms: a morphing wing and a geometric inverter.A continuous material distribution problem following the SIMP approach was formulated to find the optimal layout throughout a design domain.For our purpose, a hyperelastic Neo-hookean material was chosen.
We considered two TO problems.The first one consists in the minimization of the volume fraction of a NACA 0033 airfoil, which must be conceived as capable of reaching a fixed displacement at its trailing edge.The second TO problem aims at maximizing the geometric gain in a structure with given volume fraction, in order to design a geometric inverter starting from a rectangular initial domain.The stress distribution was considered as the main constraint in all the problems presented in this paper.In particular, a relaxation and aggregation method was adopted to treat this constraint optimally.
The problem of morphing wing design was solved using two different models of the relation between Young Modulus E and density ρ of each element, namely classical SIMP and modified SIMP (see Eqs. ( 2)-( 3)).The objective is always to minimize the volume fraction V , under constraints of performance (in terms of trailing edge displacement) and of feasibility (in terms of stress).Results showed a substantial gain in terms of V when we adopted the modified SIMP rather than classical SIMP: we pass from a volume fraction of 51.64% to 19.62%.The second structure is also easier to manufacture using a 3-D printing technique or additive manufacturing, which are the natural outputs of Topology Optimization.
Modified SIMP approach was applied to parameterize the domain in the geometric inverter problem.The first results show a ratio of 318% between output and input displacement (respectively U out and U in in Fig. 2).Several analyses were performed, by modifying the stress limits or the applied load.One can observe that increasing the stress constraint does not imply a significant gain in terms of structure performances.By analyzing the effects of the applied load, the geometric gain increases as the force increases, but the relation is not linear.
Looking at the number of iterations, by reducing the allowed stress, convergence results to be slower; on the other hand, a general dependence between number of iterations and applied load cannot be proved.
Globally, the present work demonstrated how a light structure can be easily designed through TO even with hyperelastic material responses and stress constraints.Nonlinear stress constraints were treated trough aggregation based on the lower bound Kreisselmeier-Steinhauser function: effects of different aggregation functions could be further investigated.
The position of the input force was considered as an input of each optimization problem presented in this paper, but could be treated as a new design variable as well.
Moreover, the two optimization problems could be extended to 3D analysis, allowing also the introduction of more realistic Boundary Conditions.
Future works will focus on the application of stressbased TO on the design of a morphing winglet subject to very high deformations.
Looking forward, further attempts with different optimization problems could prove beneficial to literature on the subject.

-
Central Bulb: It absorbs the applied force and transmits it to the rest of the structure.-Support Net: It guarantees an efficient redistribution of stresses all over the structure: the ramifications create the necessary moments to equilibrate the whole mechanism.In fact, as the available volume fraction constrained, the resulting positioning produces the best compromise.-Angular chain: This connects the Support Net with the clamped point, guaranteeing the force equilibrium of the structure.-Terminal chain: It make the morphing nature of the whole mechanism possible, generating the displacement in output.

Table 1 .
Default parameters used in the Airfoil optimization problem.

Table 2 .
Default parameters used in the Inverter optimization problem (symmetrical problem).

Table 4 .
(24)rical results of the Inverter optimization problem(24)with variations in F app .