Sequential topology and shape optimization framework to design compliant mechanisms with boundary stress constraints

We present a sequential topology and shape optimization framework to design compliant mechanisms with boundary stress constraints. In our approach, a density-based topology optimization method is used to generate the configuration of the mechanisms. Afterwards, a node-based shape optimization is invoked to obtain an exact boundary representation. A specialized, optimality criteria-based design update is formulated for the shape optimization. To avoid impractical hinges with point connections, stress constraints are imposed. The stress constraints are imposed using two strategies: Local stress constraints on the nodes of the boundary or global P-norm stress constraints in the domain. Further, an adaptive shape refinement strategy is adopted to increase the design space of shape optimization and to capture the fine-scale details of the geometry. Finally, numerical experiments are presented, showing that the proposed approach can be effectively applied to the design of compliant mechanisms with stress constraints.


Introduction
Compliant mechanisms are used to transmit force by elastic deformation. The design of compliant mechanisms poses a unique challenge, it requires the structure to be stiff enough to support the applied loads and also to be flexible enough to meet the kinematic requirements. Among several methods to design such a mechanism, topology optimization is the natural choice. However, within the context of topology optimization, there is no universally accepted formulation to design compliant mechanisms (Deepak et al. 2008;Zhu et al. 2020). One of the key reasons for so many formulations is to avoid narrow hinges (point flexures) in the resulting mechanism. The stresses in these hinges often exceed the material stress limit.
A mechanics-based approach to avoid narrow hinges is to impose stress constraints in the design problem. The use of stress constraints in the design of compliant mechanisms has been employed by De Leon et al. (2015). Here, a global p-norm of the von-Mises stress is used as constraint. Although narrow hinges were eliminated in the design by employing stress constraints, the authors expose some challenges. The imposing of stress constraints does not yield mesh-independent design. The method was then extended to geometric and material nonlinearity in De Leon et al. (2020). The use of nonlinear models helps to represent the real world situation more aptly. However, in both the works, the hinge region contained gray material (not fully solid cells) and lacked the exact representation of the geometry. de Assis Pereira and Cardoso (2018) presented two strategies of imposing the limit on stresses: local and global stress constraints. In their work, it was concluded that the best strategy to obtain hinge-free mechanisms is the combination of proper size control by a filter radius, a refined mesh and consideration of stress constraints. Even here, the results have gray elements around the hinges and lack the exact geometric representation of the structure. Lopes and Novotny (2016) considered a compliant mechanism design under stress constraints based on topology derivatives and a level set domain representation. The proposed approach was able to control the undesirable narrow hinges in the mechanism. In the context of level set based methods, Chu et al. (2018) presented a multi-material optimization of compliant mechanisms with stress constraints. The results show that multi-material structures without undesirable hinges are obtained and the stress constraints in different materials are simultaneously satisfied. In Emmendoerfer et al. (2020), local stress constraints were considered in the context of level set method for the design of compliant mechanisms. Here, an augmented Lagrangian approach and level set implicit boundaries were used to handle local stress constraints. The authors were able to obtain compliant mechanisms without point flexures.
A robust formulation was applied to design compliant mechanism with stress constraints by da Silva et al. (2019). This approach helps in tackling manufacturing uncertainties. Here, a stress failure criterion was considered in each projected field, in order to ensure that compliant mechanisms satisfy the stress failure criterion even in the presence of uniform manufacturing variations. It was found that the traditional deterministic approach was non-robust with respect to uniform boundary variations when compared to the robust formulation approach. The robust formulation was further extended to nonlinear elasticity along with a path-generation formulation by da . The nonlinear analysisbased optimized structure was able to provide solutions with good performance in situations of large displacements. The path-generating formulation was able to provide solutions that follow the prescribed control points, including stress robustness.
All the aforementioned works employ density-based derivatives to arrive at the desired design. Recently, there is a growing number of publications that combine shape and topology optimization techniques to obtain a design with crisp boundaries. Especially in the context of stress constraints, the maximum stresses are on the boundary of the structure. Hence, it is very important to use shape sensitivity information to control the boundary variations.
There are different approaches to combine shape and topology optimization. Eschenauer et al. (1994) introduced a bubble method, in which design updates are mainly driven by shape optimization. When the shape updates do not result in further improvement of the design, a topological change in form of a spherical hole (bubble) is introduced. This increases the design space and the shape optimization process continues. Christiansen et al. (2014) and Lian et al. (2017) applied Deformable Simplicial Complex (DSC)-based topology and shape optimization. The method represents the surface explicitly and discretizes the domain into a simplicial complex which adapts both structural shape and topology. Riehl and Steinmann (2015) used a staggered approach for optimization. The domain variation is specified by modification of boundary nodal points. Once the minimum topological sensitivity is no longer encountered at the design boundary, a hole in the domain is generated using the topological sensitivity and the procedure continues for the newly established design boundary. Stankiewicz et al. (2021) coupled topology and shape optimization into a single optimization problem. They exploit the embedding domain discretization (EDD) method to couple the shape and topology update. The structured embedding domain is assigned a pseudo-density field which is updated based on topology sensitivity. The boundary of the embedded body serves as a variable shape, which is updated by the shape sensitivity. In their work, it can be noticed that the initial design improvements were mainly driven by the topology updates, and in the later stage, the shape optimization drives the improvements. This motivated the present work to use the topology and shape updates sequentially.
Sequential use of topology and shape optimization is favorable due to their complementary nature. Topology optimization is robust in synthesizing the design from a simple initial configuration, e.g., a rectangle in the twodimensional case. However, the result of topology optimization has typically rough boundaries and does not represent the exact geometry of the structure. In contrast, shape optimization gives an exact geometry of the structure. However, it is not robust enough to arrive at complex designs starting from a simple initial configuration. In our framework, we sequentially use topology optimization followed by shape optimization to have the best of both approaches. This sequential framework is used to design compliant mechanisms with stress constraints.
In this article, we tackle the challenge to obtain a compliant mechanism with an explicit geometric description starting from a simple initial configuration (e.g., a rectangle in two dimensions). Topology optimization is employed to generate a rough mechanism layout. Next, a node-based shape optimization is performed using the embedding domain method. In the shape optimization, we impose stress constraints to obtain a durable design. To impose the local stress constraints, the stresses must be evaluated on the boundary. We present a novel procedure to evaluate the stress on the boundary in the embedding domain setting. To handle the stress constraints in shape optimization, an augmented Lagrange sub-problem is defined. To solve the augmented Lagrange sub-problem, we have novel adaptation of optimality criteria method for shape optimization.
The paper is organized as follows. In Sect. 2, the sequential framework to design complaint mechanisms is presented. Here, we explain the topology and shape optimization problems used to design the mechanisms. In Sect. 3, the results are presented. To illustrate the method, a compliant inverter and gripper benchmark problems are solved. The paper is concluded with some remarks in Sect. 4.

Method
A possible goal of the compliant mechanism design problem is to maximize the displacement at the output port, for a given input force, subjected to a restriction on the amount of material, while keeping the stresses under a certain limit. Such a design problem is shown in Fig. 1. The design domain is . A Dirichlet boundary condition is prescribed on D . The input port in is subjected to an input force f in . The output port is represented by out where the displacement should be maximized. The actuator stiffness at the input port is represented by the spring k in , and the workpiece stiffness is represented by the spring k out . The design framework involves three steps: (a) Topology optimization; (b) Shape generation; (c) Shape optimization. In this section, these three phases are presented.

Topology optimization
Topology optimization distributes material in a fixed design domain to meet the design criteria. A density-based approach for topology optimization is employed in this work. In the density-based approach, the design domain is discretized into finite elements, and each element is associated with a relative density varying from 0 (represents void) to 1 (represents solid). In our implementation, we use the modified SIMP approach (Sigmund 2007).
The problem definition for topology optimization of the compliant mechanism is given as follows: where is the variable pseudo-density (design variable), u is the displacement field obtained by solving the linear elastic boundary value problem, L is a unit vector, V is the restriction on the allowed volume of the material, e is the density corresponding to the eth finite element, and N e is the number of finite elements (design variables).
The problem in Eq. 1 is solved using the optimality criteria optimization algorithm, as outlined in Bendsoe and Sigmund (2013). Further, to have a minimum feature size and to avoid mesh dependency, we perform sensitivity filtering as explained in Sigmund (2001).

Shape generation
After topology optimization, the boundary of the optimized configuration has to be extracted. The results of the topology optimization are in a pixilated form. To obtain a geometry from the pixilated result of topology optimization, we use the alpha shape generation algorithm (Akkiraju et al. 1995). In our implementation, we resort to the alpha shape generation implementation from the CGAL library (The CGAL Project 2021; Da 2021).
A typical pixilated configuration after the topology optimization is shown in Fig. 2a. To generate a shape around the solid cells, we first need to extract the centers of all the cells above a certain density threshold as shown in Fig. 2b. Next we generate an alpha shape enclosing these points as seen in Fig. 2c. As we can see, the obtained shape is not smooth. It has jagged edges due to the pixilated nature of the topology optimization results. At this point, we compute the nodal averaging vector, which for each node, contains a displacement information necessary to reach an average position of their respective adjacent nodes. We then move the shape by a step length of 0.5 along these nodal averaging vector to obtain a smooth surface as shown in Fig. 2d.

Shape optimization
We incorporate the shape optimization method using an embedding domain (EDD) technique, introduced by Riehl and Steinmann (2017). In this approach, the physical domain boundary represents the variable shape. The shape is discretized into segments. The nodes of these segments are the design variables for the shape optimization.
In the EDD method, the variable shape of the structure is embedded within a uniform finite element background Fig. 1 Compliant mechanism design problem mesh, which is then used for the solution of the physical state problem. An adaptive mesh refinement-based boundary tracking procedure is adopted to accurately represent the structure. After the boundary tracking, the cells are categorized as interior, exterior, and boundary cells based on their relative position with respect to the physical domain boundary. Further, while solving the BVP, the exterior cells are excluded from the computation. The interior cells are treated as traditional finite element cells. For the boundary cells, we use a Gauss point oversampling integration rule, where each Gauss point is checked whether it is inside or outside the shape. If it is inside the shape, it has a standard contribution to the system matrix. If the point is outside, only a weak contribution is added to the system. The Dirichlet boundary conditions are imposed by a penalty function. The Neumann boundary conditions are imposed by performing line integration over the shape. A detailed explanation and implementation aspects of EDD are presented in Riehl and Steinmann (2017) and Stankiewicz et al. (2021).

Shape optimization strategy
Based on the fact that a local stress constraint is introduced to the design problem of compliant mechanisms, a large number of constraints must be considered in the optimization problem. To deal with it, we employ the Augmented Lagrangian formulation (Emmendoerfer and Fancello 2014; Da Silva and Cardoso 2017; de Assis Pereira and Cardoso 2018). To illustrate the Augmented Lagrangian formulation, consider a generic shape optimization problem: where f ( ) is the objective function, g i ( ) is the ith inequality constraint, and N g is the number of inequality constraints. The Augmented Lagrangian sub-problem for the optimization problem in Eq. 2 is defined as follows: where i are the Lagrange multipliers associated with the g i inequality constraint, r is a penalty parameter and the Macaulay brackets operator, ⟨⋆⟩ is defined as max(⋆, 0).
After solving the jth sub-problem, the Lagrange multipliers i and the penalty r are updated according to Steps involved in shape generation: a Density distribution after topology optimization. b Centers of solid cells for alpha shape generation. c Alpha shape enclosing the given points. d Alpha shape after smoothing by nodal averaging The shape sensitivity of the augmented Lagrangian, L is obtained on the continuum level and transformed into boundary integral form (Choi and Kim 2004). To avoid jagged boundaries and mesh dependency, we have adapted the traction method by Azegami and Takeuchi (2006). We use the traction method adapted for the embedded body (Stankiewicz et al. 2021) and geometric smoothing. Finally, to represent the shape accurately, we use the adaptive shape refinement strategy proposed by Stankiewicz et al. (2021). The adaptive shape refinement increases the design space of shape optimization and helps capture fine-scale details of the geometry.

Shape update by Optimality Criteria Method
The sub-problem in Eq. 3 is solved using the Optimality Criteria Method (OCM). For the fundamentals of the OCM, one can refer to Rozvany (1992Rozvany ( , 2012. Here, we present a brief description of our adaptation for shape optimization.
We start with defining the Lagrange function for the optimization sub-problem in Eq. 3.
where is the Lagrange multiplier for the volume constraint. The optimality condition with respect to variations of the domain is The update scheme is formulated such that the above conditions are satisfied. The design variables, i.e., vertices x on the boundary, are updated along the descent direction of the Lagrange function, i.e., Eq. 7. To have a stable smooth design update, a move limit on the design update is introduced through the traction method. This controls how much a vertex is allowed to update in each step. The move limit is introduced as a constraint in the traction method, refer to the appendix Sect. B for detail. However, is still unknown in Eq. 7. Hence, the value of is computed using Eqs. 8 and 9. Using these two conditions, a line search is performed along the search direction provided by Eq. 7. The update scheme is summarized in Algorithm 1. To perform the line search, the BVP is not (5) r j+1 = r j , > 1.
Stationary condition , Primal Feasibility condition , Complementary condition , Dual Feasibility condition , solved and the gradients are not computed again, which is a major advantage of this method.

Stress evaluation on the boundary
To impose the local stress constraints on the vertices of the shape boundary, we have to evaluate the stress values at these vertices. Due to EDD solution of the BVP, the stress values should not be interpolated from the shape function of the BVP. This is, because the accuracy of interpolation depends on the relative position of the boundary cell with respect to the shape. To illustrate the dependency, we consider two scenarios presented in Fig. 3a, b. In first case, the cell is only partially inside the shape, as shown in Fig. 3a.
The solution field in this cell will have larger error due to its weak stiffness according to the EDD method. Hence, if we evaluate the stress based on shape function interpolation, then the stress values will not be accurate. While, in the second case, a large area of the cell is within the shape as shown in Fig. 3b, the evaluation of the stress based on shape function interpolation will have better, but still not satisfactory accuracy.
A robust way to evaluate the stress on the shape is to use the L 2 projection (Zienkiewicz and Zhu 1987) of the stresses on the shape. The right-hand side of the projection system is formed by considering the contribution of all the domain cells cut by the shape segments. As illustrated in Fig. 3c, each shape segment is further split into subsegments cut by the edges of domain cells. The nodes of the subsegments are indicated by squares in Fig. 3c. On each of these segments, two Gauss points are considered and stress values are evaluated by shape function interpolation. The contributions are assembled at the respective vertices of the shape segments.
The stress values form the L 2 projection are then used to impose the local stress constraints on the boundary.

Local stress constraint design problem
We are now ready to address compliant mechanisms with local stress constraints on the boundary, which are defined as follows: where VM p is the von-Mises stress at the p th design node, a is the prescribed allowable stress, and N p is the number of points.
The Augmented Lagrangian functional of the problem in Eq. 11 without the volume constraint is The computation of the sensitivity of the above expression is presented in the appendix in Sect. A.1.

Compliant mechanism with global stress constraint
The compliant mechanism design problem with global P-norm stress constraint in the domain can be defined as follows: As before we define the Augmented Lagrangian formulation of the problem in Eq. 13 without the volume constraint, The computation of the sensitivity expression is presented in the appendix in Sect. A.2.

Summary of the novelties of the method
We have presented all the components involved in a sequential topology and shape optimization of a compliant mechanism with stress constraints. The novelty of this framework is a sequential approach to exploit the best features of both topology and shape optimization to obtain the final optimized design. Topology optimization is used to obtain an initial design. The alpha shape generation method is then utilized to obtain the boundary of the geometry. Finally, shape optimization using the embedding domain method is employed to fine tune the design.
To impose the stress constraints in shape optimization, the stresses must be evaluated on the boundary. We have presented a novel procedure to evaluate the stresses on the boundary while using an embedded domain method. To handle a large number of stress constraints in shape optimization, an augmented Lagrange sub-problem is defined.
A novel adaptation of OCM to shape optimization is presented to solve the sub-problem. Traditionally, OCM is not suitable for shape optimization since it requires the volume constraint to be satisfied. In most traditional shape optimization problems, we start with an initial configuration such that the volume constraint is not satisfied, and hence, OCM was unsuitable. However, in our framework, the topology optimization generates an initial design with volume constraint satisfied. This enables us to adapt the simple and efficient OCM to shape optimization.

Fig. 3 Stress evaluation on the boundary
Page 7 of 16 180

A note on implementation
Our program is built on the open-source finite element library deal.II, v. 9.2 (Arndt et al. 2020). We have implemented topology optimization based on Sigmund (2001Sigmund ( , 2007. At the end of topology optimization, a point cloud of all the solid cells is extracted. This point cloud is the input for the alpha shape generation algorithm. We use the alpha shape generation algorithm implementation from the CGAL library (The CGAL Project 2021; Da 2021), which is integrated into our code. The result of the alpha shape generation is a co-dimensional mesh (surface) that will be an input for shape optimization. The shape optimization using the embedding domain method is based on Riehl and Steinmann (2017) and Stankiewicz et al. (2021), and it is implemented using deal.II in the same code. All the above-mentioned steps require a tailor-made implementation in our code.

Results
We In both examples, we start with topology optimization to maximize the displacement at the output port. Starting from the configuration obtained from topology optimization, we subsequently perform shape optimization with stress constraints. The shape optimization is performed for different values of allowable stresses, and the results are presented.
The following input parameters are used throughout the result section: L = 100 in Figs. 4 and 15, Young's Modulus E = 1 , Poisson's ratio of = 0.3 , applied load of f in = 1 , and input stiffness of k in = 1 . The output spring stiffness is considered as follows: (a) k out = 0.001 , for the inverter problem, and (b) k out = 0.005 , for the gripper problem. Both mechanisms have a volume constraint of V = 0.25V 0 .

Force inverter
The first example is a force inverter. In this mechanism, the output port is expected to move in the opposite direction to the force at the input port. The problem setup is shown in Fig. 4. For topology optimization, we have used a mesh of 100 × 50 elements and a filter radius of 3.2 to design the compliant mechanism. In the present sequential approach, we can use a coarse mesh in the topology optimization phase, as it is used to only obtain a configuration of the mechanism. Further, we terminate the topology optimization once the percent of gray cells falls below 10%.

Topology optimization
The configuration of the mechanism after the topology optimization is shown in Fig. 5. To generate the boundary of the optimized configuration, the alpha shape generation algorithm is invoked. The alpha value in the alpha shape generation is set to the dimension of the element used in topology optimization. The generated shape is shown in Fig. 6.

Local stress constraint
We performed shape optimization of the force inverter with local stress constraints on the vertices of the shape   boundary. To evaluate the influence of the stress constraints, three different allowable stresses were prescribed. The resulting geometry of the optimized structure for different allowable stresses is presented in Fig. 7. It is observed that as the allowable stress increases, the thickness of the hinge decreases. Further, it is evident from the plot in, Fig. 8 that there is an inverse relation between output displacements and allowable von-Mises stresses.
In Fig. 9, the distribution of the shape vertices can be seen. We can clearly see the concentrated distribution of nodes at the regions of high curvature of the shape. The use of the adaptive shape refinement scheme allows us to represent the shape precisely.
Finally, the distribution of stresses in the structure is shown in Fig. 10. A homogeneous distribution of stresses is observed, and they are within the allowable limits.

Global P-norm stress constraint
We have performed shape optimization of the force inverter with a global P-norm stress constraint in the domain. In the present example, P = 15 is chosen. Similar to the local stress constraint case, to evaluate the influence of the stress constraint, three different allowable stresses are prescribed. The resulting geometry of the optimized structure for different allowable stress is shown in Fig. 11. Analogous to the local stress constrained structure, it is observed that as the allowable stress increases, the thickness of the hinge decreases. Additionally, it is evident from the plot in, Fig. 12, that as the stress constraint is relaxed, the objective is further minimized.
In Fig. 13, the comparison of the shape designed with different allowable stresses is depicted. Here, it can be noticed that the higher the allowable stress, the farther the hinge is moving diagonally away from the input force. Here, the shape optimization is not only modifying the thickness of the hinges but also moving the location of the hinge to achieve the desired performance.
Finally, the distribution of stress in the structure is shown in Fig. 14. It is evident that the stress is within the allowable limits.
At this point, we compare the performance of local stress constraints vs a global stress constraint. A summary of the comparison is presented in Table 1. It can be observed that the mechanisms with global stress constraints have slightly higher stresses when compared to the mechanisms with local stress constraints. In both cases, it is clear that as the allowable stresses increase, the performance of the structure is better.

Gripper
The second example is a gripper mechanism. The intent of this mechanism is to transform an input force in horizontal direction to an output force in vertical direction. The problem setup is shown in Fig. 15. We use a mesh of 100 × 50 elements and the filter radius of 3.2 to design the gripper mechanism. We terminate the topology optimization once the percent of gray cells falls below 10%.

Topology optimization
The configuration of the gripper after the topology optimization phase is depicted in Fig. 16. The alpha shape generation algorithm is then invoked to generate the shape of the gripper. The generated shape is shown in Fig. 17.

Local stress constraint
Similar to the previous example, we now perform shape optimization of the gripper with local stress constraint on the boundary nodes. To study the influence of the stress constraint, three different allowable stresses are prescribed. The geometry of the optimized structure for different allowable stress is shown in Fig. 18. The inverse relation on the hinge thickness and the value of allowable stress can be observed in the figure. Also in the plot in Fig. 19, a compromise can be observed in the prescribed allowable stress and the objective.
The stress distribution in the gripper is shown in Fig. 20. An even distribution of stress is observed and the stresses in the structure are within the allowable limits.

Global P-norm stress constraint
At this point, shape optimization of the gripper with a global P-norm stress constraint in the domain is performed. In the present example, P = 15 is chosen. Similar to the local   stress constraint case, we prescribe three different allowable stresses to evaluate the influence of the stress constraint. The optimized structure for different allowable stress is captured in Figs. 21 and 22. Further, it is evident from the plot in Fig. 12 that as the stress constraint is relaxed, the objective is further minimized. Finally, a comparison of the performance of local stress constraints vs a global stress constraint is performed, and a summary is presented in Table 2. In this case, it can be observed that the mechanisms with global stress constraints have violated the allowable stress by a significant margin.

Conclusions and outlook
A sequential topology and shape optimization framework is developed to design compliant mechanisms. Topology optimization is robust in generating an optimal configuration and, however, does not give an exact geometry of the structure. Shape optimization proves to be effective at fine tuning the design and giving the exact geometry of the structure; however, it is not robust enough to generate optimal configurations. In our approach, we incorporate these methods to get the best of both worlds. We start with topology optimization to generate the configuration and continue with shape optimization to generate fine-scale detailing and obtain the exact geometry. In the process, we developed a novel optimality criteria-based design update for shape optimization with volume constraint. Further, the use of an adaptive shape refinement strategy in shape optimization helps us expand the design domain and capture the fine-scale detail of the geometry.
Finally, we presented two benchmark examples of compliant mechanisms. In these examples, we designed mechanisms to maximize the displacement at the output port with a restriction of the material use and stress constraints. The stress constraints were imposed using two strategies: (a) The methodology presented in this work is restricted to linear elasticity. Problems like the displacement inverter might require consideration of large displacements, i.e., geometric nonlinearity. Moreover, the availability of exact boundary representation in the shape optimization step offers great potential for development of geometric or manufacturing constraints, which are otherwise difficult to implement in topology optimization, e.g., curvature radius constraints for milling.
The presented framework can in principle be extended to three dimensions. Topology optimization is well established in three dimensions, and its major drawback is its computational cost. In our framework, fine geometrical details are captured by shape optimization; thus, a coarser grid for topology optimization can be employed. Moreover, shape generation by the alpha shape generation algorithm is in principle available for generating bounding surfaces of a three-dimensional given point cloud. Alternatively, there are also advanced tools in computer graphics to generate surfaces given a point cloud from topology optimization. However, shape optimization using the embedding domain method in three dimensions is most challenging. Therein, the major obstacles are the mesh tracking and sub-triangulation of faces intersecting with embedding cells, which come with a high computational cost. Taking together, the implementation in three dimensions is in principle feasible but requires further investigation.

Continuum approach to shape sensitivity
Shape sensitivity analysis using the continuum approach is performed based on Choi and Kim (2004). The material derivative concept is applied, which is not explained her, but interested readers can refer to Arora (1993); Choi and Kim (2004). We first get the sensitivity expression for a general response functional based on the presentation  in Stankiewicz et al. (2021). Consider a general response functional We apply the adjoint method and augment the above functional by a weak form of equilibrium equation. The variation of the state field in the weak form is replaced by the adjoint variable: where Using variational calculus, it was proven that W a = 0 (Arora and Cardoso 1992). This leads to We then apply the boundary integral method presented by Arora (1993) to represent the material derivative of augmented functional: The partial derivatives of the integrands Ḡ , ḡ, and h are as follows: Introducing A.15 and A.17 into A.19 and grouping of terms gives us: It is required that the variation of the augmented functional should vanish according to variational principle for design sensitivity analysis (Arora and Cardoso 1992). This means . (A.20) .
. the expressions in the two curly braces should vanish. The expression in the first curly braces in Eq. A.21 represents the weak form of the equilibrium for the primary structure, with the explicit derivative of the adjoint state field [u a ] � serving as the test function. So this term naturally vanishes. Further, the expressions in the second curly braces in Eq. A.21 are used to define the adjoint problem such that even this term vanishes. In the adjoint problem, we utilize the symmetry of the Cauchy stress tensor and the major symmetry of the elasticity tensor in which u ′ takes the role of admissible test function. Leading to the following definition of the adjoint problem, With this, the final sensitivity expression is

Complaint mechanism with local stress constraint
The augmented functional for compliant mechanism with local stress constraint from Eq. 12 is To obtain the adjoint problem, we substitute Eq. A.24 in Eq: A.22, which leads to The sensitivity expression is obtained by using Eq. A.24 in Eq: A.23, .

Complaint mechanism with global stress constraint
The augmented functional for compliant mechanism with global P-norm stress constraint from Eq. 14 is To simplify the notation, we define the following placeholders: Now the adjoint problem is obtained by substituting Eq. A.27 in Eq: A.22, which leads to .

Traction method with move limit
We use the traction method presented by Stankiewicz et al. (2021), which is adapted for embedding domain discretization. Briefly explained, the traction method is based upon solving an auxiliary BVP, in which the raw shape sensitivity is employed in form of external, nodal traction forces. The solution of such BVP yields a smooth, mesh-independent design update vector. To the traction method in the aforementioned article, we have added an additional energy term to control the move limits. The new energy form with the penalty for move limit is given by where ext is the energy of the external forces, i.e., the raw sensitivities, P is the penalty functional, which is responsible for the bounding box constraints, P ml is the new penalty to control the move limits in each iteration, and is the scaling factor which ensures that the solution of the traction method BVP yields an exact design update vector. The definition of the move limits penalty functional is where n c is the number of points subjected to move limit constraint, p g c is the penalty function of which values are dependent on the constraint violation measure g c which are defined in below equations: where ū is the allowed shape change. To enable a smooth design update, a move limit of 0.01 was chosen for all the simulation presented in result section. (A.31) . (B.32) t = int + ext + P + P ml → Min,