Stress-based shape and topology optimization with the level set method

This paper proposes a level set method to solve minimum stress and stress-constrained shape and topology optimization problems. The method solves a sub-optimization problem every iteration to obtain optimal boundary velocities. A p- norm stress functional is used to aggregate stresses in a single constraint. The shape sensitivity function is derived and a computational procedure based on a least squares interpolation approach is devised in order to compute sensitivities at the boundaries. Adaptive constraint scaling is used to enforce exact control of stress limits. Numerical results show that the method is able to solve the problem e ﬃ ciently for single and multiple load cases obtaining solutions with smooth boundaries.


Intro
Stress optimization is one of the key factors for structural design in a wide range of engineering problems. Its importance relies in obtaining mechanical structures that are free of or present less stress concentrations. Topology optimization has shown to be a powerful tool for structural conception. If stress is not taken into account in the concept of the structural component, the design is susceptible to postprocessing or rework, leading to unexpected costs. However, considering stress within the context of topology optimization has been an arduous task for our research field since the pioneering paper by Duysinx and Bendsøe [1].
Notably, Le et al. [2] presented a practical solution for stress-constrained design in the framework of density-based topology optimization [3]. The authors used a set of numerical ingredients to circumvent three major issues namely the singularity problem, the local nature of the stress and the highly nonlinear stress behavior. The first refers to singular optima and the difficulties that nonlinear programming algorithms have to identify degenerate subspaces where the globally optimal stress design can be located. In practice, very low density elements can present high stress values, leading to the non-removal of such elements. Stress relaxation has been applied to ensure that the element stresses and volume tend to zero simultaneously when the latter approaches zero [4,5,6]. The second issue occurs because stress is a local quantity, which implies that stress levels must be controlled in every material point of the structure. This leads to a large number of constraints. Global measures of the maximum stress have been used to estimate and control the stresses throughout the structure. Such measures can be based on the p-norm or Kreisselmeier-Steinhauser (K-S) functions [2,7,8,9,10,11]. The third related issue is that stress levels are highly sensitive to design changes. This can be specially critical in stress concentration regions, such as sharp and re-entrant corners.
The breakthrough in Le et al. [2] has been the normalization of a p-norm stress measure to approximate the peak stress in the structure, applied together with stress relaxation and consistent density filtering in a solid isotropic material with penalization (SIMP) framework. A key factor to the success of this method is that it uses a consistent density filtering versus sensitivity filtering, the latter mentioned to not perform well because of the nonlinearity of the stresses. Since then, stress-constrained topology optimization has attracted a lot of attention from researchers, including different formulations and models such as a damage approach [12], fluid-structure interaction [13] and thermo-mechanical coupling [14].
In the level set optimization framework, stress was first considered in the papers by van Miegroet and Duysinx [15] and Allaire and Jouve [16]. The first carried out shape optimization to alleviate stress concentrations in 2D fillets. The later presented stress shape derivatives and minimum stress topology solutions. One advantage of the level set method (LSM) is that the structural boundaries are well defined throughout the optimization [17].
An extensive range of techniques have been proposed throughout the last years to incorporate stress-based design into level set topology optimization methods. Most of the papers propose new penalty functions or global stress measures [18,19,20,21]. Verbart et al. [22] employed in the LSM the same set of techniques previously used in density-based methods, such as stress relaxation, p-norm aggregation function, a SIMP model and adaptive constraint scaling. Xia et al. [23,24] enhanced the performance of piezoresistive sensors by tailoring the performance of stress in regions of interest via a level set based method. Topological derivatives for stress-based functionals were also developed [25,26] and recently applied to topology design of compliant mechanisms [27]. Wang and Li [28] formulated stress-constrained topology optimization with a novel shape equilibrium problem of active constraints, removing the intrinsic nondifferentiability introduced by local stress constraints. The same idea was applied later into a stress isolation problem [29]. Activation and deactivation of constraints were also applied by Emmendoerfer and Fancello in [30] and [31], together with a few regularization techniques. The latter paper used a reaction-diffusion equation to guide the design optimization sequence in order to eliminate the level set reinitialization steps. It was shown to improve stability and convergence and new holes inside the domain were created, although the stressconstrained results did not confirm the low dependency of the initial level set function.
Most of the works mentioned above apply the Finite Element Method (FEM) or the extended FEM (XFEM) and rich discussions arise regarding the accuracy of stress computations in discretized domains. In the same scope, Cai et al. [32] and Cai and Zhang [33] combined the Finite Cell Method (FCM) with spline functions to model any type of geometry and enhance stress accuracy. Guo et al. [34] used the XFEM to develop stress related optimization with multi-phase materials. Although considerable progress was achieved, most of the works relating stress and level set topology optimization are not able to enforce stress constraints without the aid of extra numerical techniques, e.g. perimeter regularization or a compliance term added to the objective. Other papers present oscillatory boundaries and/or very thin structural members in their solutions, most of them not being able to actually remove the re-entrant corner in the benchmark L-bracket example. More recently, Poljnar et al. [35] presented good topology solutions for stress minimization and proposed a novel objective based on a penalized stress-deviation measure and approximated the level set description by the XFEM. The improvements in level set based methods are highly significant to problems in which boundary location is important, such as in stress and many others, like multiple-type boundary problems [36,37].
This paper presents the level set topology optimization method for stress minimization and stress constrained structural design problems. The classical finite element method with ersatz material is used. Within the context of stress-based level-set based methods, the contributions are the following: • The use of a p-norm function to aggregate stress constraints. A p-norm stress functional was previously applied to a level set topology optimization in [26] and in [22]. However, the first used a weighted compliance to counter balance the nonlinearities of the p-norm. The latter used a level set function but with all densitybased optimization numerical ingredients. Herein, we apply a pure p-norm stress functional and their shape derivatives.
• A least-square interpolation is used to obtain a more accurate estimation of the stresses used in the calculation of the sensitivities.
• The velocity field is constructed by a Lagrangian function with shape sensitivities of the objective and constraint functions, in which an optimization problem is solved at every iteration to determine the optimal multipliers, 2 reducing the optimization variables to a set of one multiplier per objective and constraint.
• The adaptive scaling from [2] is used to ensure exact stress limits in stress-constrained problems. The work in [22] was the only one that previously applied this constraint scaling scheme to a level set topology optimization.
The main advantage of our method is the generality of good quality topology solutions when compared to other stress-based level set methods in a range of problems such as minimum stress, stress constraints and multiple loads and stress criteria cases. Figure 1 presents the L-bracket designs presented by several stress-based level set optimization works. The majority of the designs consider the load at the midpoint of the right-hand side edge. This makes convergence easier because it leads to a load path that already tends to create a curve around the corner. Comparing our results only to the works that put the load in the corner, our method can present more rounded shapes and be twice as efficient depending on the p-norm parameter. We postulate throughout the paper that these results are possible due to a combination of numerical ingredients not used before in the context of the level set methods.
The remainder of the paper is organized as follows. Section 2 presents the level set topology optimization method used in this work. Subsequently, Section 3 details the contributions on stress-based problems. In Section 4 numerical results and discussions are presented. Finally, Section 5 concludes the paper.

Level set topology optimization method
In the level set topology optimization, the structural boundary represented by an implicit function is updated to find a topology that minimizes an objective function while satisfying specified constraints. The conventional approach defines the structure using an implicit signed distance function: where φ is the level set function, Ω is the structural domain, Γ is the structural boundary and x is the point inside the design domain, Ω d , that contains the structure. The structural boundary is optimized as the implicit function is iteratively updated by solving the following advection equation: where t is a fictitious time domain in which the level set evolves. The velocity field of motion by mean curvature is where V n is the normal velocity and N is the unit outward normal. The discretized version of Eq. (2) with Eq. (3) is solved using an explicit forward Euler-scheme as: where k is the iteration number, j is a discrete boundary point, ∆t is the time step and V n, j is the velocity at the point j, considered as an advection velocity of the type dx/dt = V n φ/| φ|.

Level set topology optimization
The velocities required for the level set update are obtained by solving an optimization problem. A generic optimization problem can be formulated using the position of the structural boundary as the design variable: where f (Ω) is the objective function and g i is the i th inequality constraint function. The objective and constraint functions are linearized about the design variables at each k th iteration using a first-order Taylor expansion: where ∆Ω k is the update for the design domain Ω andḡ k i is the change in the i th constraint at iteration k. In the level-set description of the boundary, shape derivatives provide information about how a function changes over time with respect to a movement of the boundary point. They usually take the form of boundary integrals [40]. In this case, where s f and s g i are the shape sensitivity functions for the objective and the i th constraints. In Eq. (7) and (8) the shape sensitivities are multiplied by the time step to compute the total change in the functions after updating the implicit function using Eq. (4) with the specified time step and velocity function. Discretizing the boundary at nb points, one can rewrite: where l j is the discrete length of the boundary around the boundary point j, C f and C gi are vectors containing integral coefficients and V n is the vector of normal velocities. For a constrained problem, one can write where d is the search direction for the boundary update and α > 0 is the actual distance of the boundary movement. Expressing V n as a function of α and d, the optimization problem from Eq. (6) can be written as: It can be pointed out that the superscripts k indicate that the optimization problem is solved at every iteration. The minimum (V n,min ) and maximum (V n,max ) velocities are limited to the Courant-Friedrichs-Lewy (CFL) stability condition. In the proposed method the CFL condition enforces the maximum boundary movement to be less than the spatial grid discretization, i.e., less than one element size. For Eq. (12), the Lagrangian function can be expressed as: where λ k is the vector of Lagrange multipliers for the ni inequality constraints, µ k is the Lagrange multiplier for the equality constraint in Eq. (12) computed at iteration k. When the optimality condition is satisfied, one can assure that d is a solution of the optimization problem in Eq. (12). Solving Eq. (14) for d k and substituting in the equality constraint of Eq. (12), it is obtained: The substitution of Eq. (15) into Eq. (12) eliminates the equality constraint and the optimization problem becomes: In the present method, α is limited by the CFL condition because its direct relation with the normal velocities via Eq. (11). Additionally, it is important to point out that the set of optimal velocities is defined by a Lagrangian function with the sensitivity integrals C f and C i via the λ multipliers. The function for velocities showed to produce smooth boundaries and the possibility to address multiple constraints in [38,39]. Therefore, as long as the sensitivity field provided is smooth enough, no post-processing of velocities is needed. Also, the use of the multipliers reduces the size of the problem by the same number of existent objective and constraint functions. In the present work, the values for λ are obtained via the SQP method present in the open source library NLOPT, thus computing the function for optimal velocities. The vector of optimal velocities is then substituted in Eq. (4) to update the level-set boundaries. The following section presents the proposed methodology to handle stress optimization in the context of the level set method.

Problem formulation
The p-norm is used here as a stress aggregation function. This function approximates the maximum stress in the structure and the approximated value is closer to the actual maximum as p increases. This global stress measure can be used either as objective for stress minimization, or as an inequality constraint. The numerical results section explores these different cases. The following p-norm functional defined in the domain Ω can be written: where p is the p-norm parameter and σ vm is the von Mises stress.

Shape sensitivity function
In order to solve both stress minimization and stress-constrained problems, the shape sensitivity of the p-norm functional G(Ω) is needed. For convenience, we introduce the function thus, the derivative of the integral function from Eq. (17) can be expressed as We use the adjoint sensitivity analysis method to efficiently obtain the shape derivative J of the auxiliary function. One can define a general functional and the state equation in which a is written in terms of the elasticity tensor C and the mechanical strains ε, and, where b is the body force and F is the traction applied on the boundary Γ N ⊂ ∂Ω. Therefore, the augmented functional can be expressed as 6 in which λ is the adjoint variable. The material derivative method is applied here in order to obtain the shape derivative of the augmented functional [28,41]. The derivative of the augmented functional Λ can be expressed as where, in which V n is the normal velocity component. Substituting Eq. (26), (27) and (28) into Eq. (25) and collecting all the terms that contain λ , the weak form of the state equation is recovered, whose total sum is zero. Collecting all the terms that depend on u and letting their sum be zero, we obtain the following adjoint equation where j (u) indicates the adjoint pseudo-load. Finally, assuming that body forces are not present in our problem, i.e., b = 0 in Eq. (28), and collecting all the remaining terms, one obtains the material derivative of the augmented functional with respect to the boundary points on Γ, The derivative Λ of the augmented functional expressed in Eq. (30) is exactly equivalent to J (Ω) due to the adjoint method. Recalling that the term j (u; x) in the general functional from Eq. (20) indicates σ p vm in the present problem, the total shape sensitivity of the stress p-norm function can be computed substituting Eq. (30) in the term J (Ω) of Eq. (19): Equation (31) indicates how the p-norm stress functional changes when a structural boundary point moves in its normal direction. The computational scheme for computing this function at the boundaries is outlined in the following section.

Computational procedure
The computational scheme for the structural analysis is based on the FEM on a uniform and fixed grid using bilinear quadrilateral elements. Although not necessary, the FEM and the LSM meshes are considered in this work to be exactly the same for simplicity. The level set function Φ from Eq. (1) is discretized at the grid nodes. The boundary Γ is determined using the nodal values of the level set function and discretized into a set of points that are located where the level set boundary (Φ = 0) intersects the elements edges, as illustrated in Fig. 2. The discretized level set is approximated by the finite elements. The most common method for calculating the stiffness of the cut boundary elements is to weight the element stiffness by an element density, here defined as the area of the structure within a cut element divided by the total element area. This is also known as the ersatz material approach. The formulation of the element stiffness is where K e is the element stiffness matrix, ρ e is the element density and K 0 is the stiffness of an uncut element. After assembling the global stiffness matrix, Eq. (21) can be solved for the displacement field u.
Equation (31) is derived in the continuum. In order to apply it in the proposed level-set methodology, the sensitivities must be computed in the discretized boundary points. This is carried out via a least squares interpolation scheme [42]. It starts by computing the shape sensitivities terms in the Gauss points of the finite elements, relying in the accuracy of computation in those points. The following equations express how to compute the shape sensitivities terms at a Gauss point n. The value of the functional presented in Eq. (31) is computed via the finite element discretization, as where υ i and σ i are, respectively, the volume and the von Mises stress σ vm for each element in the mesh and nel is the total number of elements. These stresses are computed at each element centroid as where σ denotes the stresses in Voigt notation and V is the matrix in a two dimensional case. The stress vector can be computed as where C is the elasticity matrix, B is the strain-displacement matrix and u the vector of displacements. The shape sensitivity from Eq. (31) is expressed in terms of the mechanical and adjoint strains. The adjoint strain is computed via the solution of the adjoint problem with pseudo-load j (u). This term can be expressed as where, and from Eq. (34), finally from Eq. (36), Substituting Eq. (38), (39) and (40) into Eq. (37) and rearranging the expression, the pseudo-load is computed as Thus, the adjoint problem integral from Eq. (29) can be expressed as where the subscript i indicates the term evaluated at the i th element and K is the global structural stiffness matrix. The pseudo-load can be assembled similarly as a body vector for each element, computed in its centroid. The final shape sensitivity of the stress p-norm function can then be computed at a Gauss point n as where the subscript n indicates the terms computed at the coordinates of the Gauss point n.
The sensitivity at a boundary point is obtained by interpolating the Gauss point sensitivities using a least squares method. All the Gauss points sensitivities inside the interpolation region around the boundary point and that belongs to elements with densities greater than 0.1 are taken into the interpolation scheme, like illustrated in Fig. 3. The sampled data are fitted to a model, which usually takes the form of a basis function. The investigations in [42] found that the following second degree polynomial provides sufficient accuracy for the sensitivity distribution in two dimensions: ζ (x, y) = L 0 + L 1 + L 2 y + L 3 xy + L 4 x 2 + L 5 y 2 , where ζ (x, y) is the sensitivity value at the Gauss point, x and y the relative distance to the boundary point in x and y directions, respectively, and L n are the unknowns computed by the least squares interpolation process. The Gauss point sensitivities and the relative distances are weighted by the element density and the inverse distance from the boundary point as: where α is the area fraction and |x p − x i | is the distance between the boundary point, x p and the sampled sensitivity value, x i . An overdetermined system of equations is then created by grouping all the sampled data (Gauss point sensitivities) in the following form: where A contains the weighted relative distances, ı contains the set of weighted Gauss points sensitivities and L n is the least squares solution. The value L 0 is the sensitivity obtained for the boundary point. As such, the sensitivity value at the boundary point is an interpolation of the sensitivity field computed inside the structure. This overcomes the difficulties in computing accurate values at the structural boundaries and the ersatz material approach can be used to provide accurate enough sensitivity for the least squares interpolation. It is important to point out that in [42], normalized velocities are averaged at the nodes, while herein no averaging is carried out. Here, we follow the investigation presented in [42] that shows that a radius of 2h, with h being the element size, consistently produces good results. Therefore, we employ the same interpolation region radius in all of our numerical experiments. The least squares method has been applied to compute the stress terms σ vm and Cε (u) · ε (λ) from Eq. (43) at each boundary point independently, before summing them up in the final shape sensitivity equation. This showed to reduce spatial oscillations in the computed velocities along the boundaries.
After computing the sensitivities and solving Eq. (16) for optimal velocities, the level set boundaries are propagated using the fast marching method [43]. In this work, the level-set function has been reinitialized at every iteration. A new set of boundary points is created in the beginning of every iteration, after boundary propagation. Convergence is reached when the following expression is satisfied for a small number τ: which evaluates the change in the objective function during the last five iterations.

Optimization procedure
Our implementation is based on an object oriented programming framework using C++ developed by Hedges et al. (2017) [44]. The optimization procedure is as follows: 1. Instantiate the finite element and the level set meshes, including material properties and optimization parameters. 2. Initialize the level set function Φ for an initial structural design, e.g., including initial holes. 3. Discretize the level set boundaries into a set of boundary points. 4. Compute area fractions. 5. Assemble and solve the finite element equations to obtain the displacements field u. 6. Solve the adjoint state Eq. (42) and compute the stress terms σ vm and Cε (u) · ε (λ) from Eq. (43) at the Gauss points. 7. Interpolate both stress terms separately to the boundary points via the least squares approach and then compute the complete sensitivities from Eq. (43). 8. Solve sub-optimization problem from Eq. (16) for optimal α k and λ k . 9. Compute the velocities at boundary points using Eq. (11).
10. Compute the level set gradients φ k j at boundary points. 11. Update level set function Φ at the grid nodes around the boundaries using Eq. (4). 12. Reinitialize level set function if necessary/desired. 13. Check convergence, e.g., the change in the objective function during the last five iterations. 14. Iterate from step 3. 10

Numerical results
This section aims to present the numerical investigation on the proposed methodology. First, a discussion on the p-norm sensitivities is presented. Then, numerical solutions for stress minimization and stress-constrained problems are presented, followed by an example with multiple stress criteria and load cases and a beam with a notch.

L-bracket
First, the benchmark L-bracket example is considered [1]. The L-bracket is modeled by a 100×100 finite element mesh with a 60×60 section removed to create the L-bracket domain, as seen in Fig. 4(a). A vertical load F = 3 units is applied. The Young's modulus of 1.0 and a Poisson's ratio of 0.3 are set as material properties. Figure 4(b) presents the von Mises stress field for such model. As devised in the previous section, the shape sensitivity is first computed at the Gauss points, obtaining then a field of domain sensitivities. Figure 5 presents field plots of the shape sensitivities. The values are raised to the power 0.25 for better visualization of their magnitude; correspondingly, they are all positive in the plot. Therefore, Fig. 5 also shows the sign of the averaged sensitivity value within each finite element. One can observe that the sensitivity fields for higher p-norm parameters present more regions with positive sign. The critical areas such as the re-entrant corner can also be clearly seen, where a zero curved isoline touches the corner. The stress values are also given in Fig.  5 for each p parameter.
The sensitivity plots show that when increasing p, more regions with positive sensitivities appear. These regions indicate a positive change in the function, i.e., an increase in the p-norm stress, if the density of that region is increased. In the level set context, when these sensitivities are interpolated at the boundaries, it indicates that an increase in stress occurs when that boundary point moves along the direction of the outwards normal vector. So, in order to decrease stress, that boundary point should move inwards the structure. This illustrates the ability of the p-norm sensitivities to reduce stress concentrations in such problems. From the sensitivity plots we can conclude that the minimum p-norm parameter to remove the re-entrant corner and reduce stress is p = 6, which is also in accordance with the literature [2]. Figure 6 presents surface plots of the domain sensitivities using p = 6. It evidences how steep are the changes in the sensitivities around the re-entrant corner of the L-bracket, which turns out to be a challenge when computing the sensitivity values at the boundaries.
We note that for p = 2, the p-norm functional is similar to the squared von Mises stress functional used by Allaire and Jouve [16]. Our solution for the same example presented in [16] but using p-norm sensitivities with p = 2 is given in Fig. 7. In this case the load is located at the midpoint of the right edge in the L-bracket, as per Allaire and Jouve [16] for comparison. The optimization is not able to remove the re-entrant corner with such sensitivities. In this manner, the global stress sum is reduced because volume is being removed, but the maximum stress is practically unaltered since the stress concentration is still present.

Stress minimization
Herein stress minimization of the L-bracket model presented in Fig. 4 is explored. In this case, the optimization problem consists in the volume-constrained stress minimization, expressed as: whereV is the limit for the volume constraint.
Starting from a fully solid design domain, Fig. 8  to circumvent an ill-posed analysis. The volume fraction limit is set as 70% for this case. The optimal solutions are hook-like structures. Since there are no initial holes, we expect only the external boundary of the L-bracket to be optimized. One can observe that as the value of p is increased, the final shape of the structure tends to be more rounded and to look more like mechanical hooks. Additionally, the final maximum stress is reduced for a higher p. Figure 9 presents snapshots of the optimization history of the level set function in this example. tural design, Fig. 10(a). This allows the boundaries of the internal holes to move and merge with each other, leading to a non-trivial and optimized structural topologies, shown in Fig. 10(b)-(d). For this example, the optimal structure was constrained to a volume fraction of 40%. Figure 11 presents the convergence history of the optimization carried out with p = 6.
Similarly to the hook-like case, as p is increased more rounded corners are obtained in the final solutions. Consequently, it can be observed that the final stress values are reduced for higher p's and the stress distributions in the final topologies are much smoother than those in the initial design. Furthermore, the manifestation of the stress optimization can be seen not only numerically but also in the shape of the boundaries. The solid regions connecting two or more structural members are larger and rounded for these solutions, increasing with higher values of p and creating structural joints that are reminiscent of the ones used in trusses.
The least squares interpolation method is used to enhance the accuracy of the sensitivities computed at any boundary point by using a set of values computed at the Gauss points. In order to demonstrate the effect of this technique, Fig. 12 presents the snapshots of the L-beam design without the least squares interpolation. In this case, the sensitivities are computed directly at the boundary points using finite element shape functions interpolation. It is observed that the solution takes the equivalent path of the ones using the proposed least squares technique. However, large oscillations in stress concentrated regions drastically reduce the convergence speed. The solution in Fig. 12 was stopped at iteration 1500 without satisfying a convergence tolerance of 0.001. The oscillations are more critical when the p value is increased. This example shows that the present method can work without the least squares interpolation, but with very slow convergence.

Stress-constrained volume minimization
For stress-constrained problems the p-norm functional itself cannot be used to enforce a constraint on the actual maximum stress σ max since, as aforementioned, the p-norm is always larger than the maximum for finite p. Herein, we employ the adaptive scaling constraint from Le et al. [2]: whereσ is the stress constraint limit in the structure and c is computed as follows: where η ∈ (0, 1] controls the variations between c k and c k−1 . The scaling uses information from the previous optimization iteration to normalize the constraint as c k σ PN so it approximates the actual maximum stress σ max as the design converges. Therefore, we formulate the stress-constrained problem as: where the objective is the structural volume. Starting from the initial design and using p = 6, the solutions obtained for different stress constraint limits are presented in Fig. 13. The initial design corresponds to iteration 0 in the figure.
It can be observed that the optimal volumes obtained are not proportional to the stress constraint limits. This can be explained by the non-monotonic behavior of the stress with respect to the volume, indicating that there might not exist an optimal volume that reduces stress to a prescribed value if material is not redistributed in the stress concentration 16 regions. In this example, for the lowest stress constraint value the optimal volume was considerably higher than the others. It can also be pointed out that for this volume minimization example, the solution obtained using the highest stress limit presented a final stress field closer to a fully-stressed design. Figure 14 presents intermediate designs for the minimum volume solution obtained using p = 6 andσ = 2.2, while Fig. 15 presents corresponding stress fields. Figure 16 presents the convergence history of the case with p = 6 andσ = 2.2. The scaling factor c k starts as 1, leading the constraint equation to be exactly the p-norm. It slowly approximates the actual maximum stress σ max based on the iteration history. It can thus be said that the scaling helps the algorithm to converge to the optimal maximum stress constraint.

Multiple stress criteria and load cases
In this section the proposed methodology is applied to solve a problem with multiple stress criteria and load cases. The T-bracket design domain presented in Fig. 17(a) is explored. The same material properties from the L-bracket design are used and F = 3 units for both load cases. This example is similar to the one presented by Le et al. [2] that allows different stress criteria to be included. The element stresses σ vm are normalized for any stress criteria based on the location of the design domain, e.g., for the T-bracket showed in Fig. 17(a) the normalization based on the coordinate x in the Cartesian system is: which allows each single load case to have multiple different stress limits in the domain, e.g.,σ 1 andσ 2 . The load cases are then grouped in one single p-norm functional as  where nlc is the total number of load cases andσ x depends on the different stress criteria in the domain Ω. Thus, applying the adaptive scaling c k from Eq. (50) we make sure that all stress criteria are satisfied with one single constraint. We then formulate the stress-constrained optimization problem as: This approach is consistent with the one in Le et al. [2]. Starting with the initial design shown in Fig. 17(b), the solution obtained for two load cases and a single stress criterion usingσ 1 =σ 2 = 1.8 and p = 6 is presented in Fig. 18. The minimal volume fraction obtained is 47.18%. For two load cases and two stress criteria, Fig. 19 presents the solution usingσ 1 = 1.4 andσ 1 = 2.2. In this case, similarly to Le et al. [2], more material is distributed to the weaker left region to accommodate the lower stress limit. The minimal volume fraction obtained is 43.11%. The normalized stress fields (σ vm /σ x ) are shown for all load cases.

Beam with a notch example
This example applies stress minimization in a beam model with the presence of a notch, similar to the one presented by Polajnar et al. [35]. The model is shown in Fig. 20(a)  For a final volume fraction constraint of 50%, the solution for compliance minimization and its correspondent von Mises stress field are presented in Fig. 21(a). Using p = 6, the solution for stress minimization and its corresponding von Mises stress field are shown in Fig. 21(b). The compliance solution presents a maximum stress of 3.3095, while the minimum stress design has σ max = 1.9673, for a reduction of approximately 41%. It is observed that stresses are better distributed in the stress-based solution and the structural members present the features described in the end of Sec. 4.1.1, such as larger rounded joints. A significant difference between the stress-driven and the compliance-driven optimal designs is the appearance in the former of wider, low stressed members at the center connecting the bottom and top members. An inspection of the sensitivity field in these members showed the sensitivities are in the same order of magnitude of those for the other members, which indicates they are equally important in reducing the high stress that occurs at the notch as other members. Therefore, these wide members draw load from the bottom member towards the top member, which helps reduce the stress at the notch.

Conclusions
This paper proposes a level set method to solve minimum stress and stress-constrained topology optimization problems. The numerical results exhibit the features of efficient stress-optimized solutions and indicate that the method is able to solve the problem effectively due to the combination of the numerical ingredients, the linearization of the optimization problem at every iteration while solving for optimal velocities, the p-norm stress aggregation, the  shape sensitivity computational procedure and the adaptive stress constraint scaling. Plots of the shape sensitivities computed in the entire structural domain indicate how the proposed stress aggregation leads to the removal of stress concentrations. Hook-like solutions were used to verify our method for shape optimization. The least square interpolation approach provides enhanced sensitivities at the boundaries and the Lagrangian function in the optimization problem efficiently leads to good quality topology solutions that compare well with solutions in the existing literature, see Fig. 1. The method is able to enforce exact stress constraints using adaptive scaling and, to the best of the authors' knowledge, it is also the first level set approach to demonstrate multiple load cases and different stress criteria in different regions of the design domain.
Volume minimization Two load cases Single stress criterion for each casē