Damage approach: A new method for topology optimization with local stress constraints

In this paper, we propose a new method for topology optimization with local stress constraints. In this method, material in which a stress constraint is violated is considered as damaged. Since damaged material will contribute less to the overall performance of the structure, the optimizer will promote a design with a minimal amount of damaged material. We tested the method on several benchmark problems, and the results show that the method is a viable alternative for conventional stress-based approaches based on constraint relaxation followed by constraint aggregation.


Introduction
Topology optimization of continuum structures has become an increasingly popular design tool since its introduction by Bendsøe and Kikuchi (1988). Although in the last decade topology optimization has found its way to industry, in most practical applications there is still a relatively big gap between the topology optimized design and the actual manufactured design. In general, a number of post-processing steps are necessary to make the design suitable for manufacturing and to meet relevant structural criteria, such as stress and buckling. Consequently, including stress constraints in the topology optimization process would reduce the necessary post-processing of the topology optimized design. However, including local stress constraints has been a major challenge because of several difficulties that arise.
First, the stress is a local state variable, which makes the problem computationally expensive. Traditionally, topology optimization has been used to solve problems of many design variables and a few responses, such as minimizing compliance under a volume constraint. These type of problems can be solved efficiently in an adjoint formulation. However, when considering stress constraints, the number of local constraints is of the same order as the number of design variables. Consequently, there is no benefit in using an adjoint formulation, and solving the problem by gradient-based optimization becomes computationally expensive.
Secondly, so-called 'singular optima' arise in stressconstrained topology optimization. Singular optima are (local) optima that cannot be reached using ordinary gradient-based optimization. Singular optima were first observed in truss topology optimization by Sved and Ginos (1968). They demonstrated on a simple three-bar truss problem under multiple loading conditions that the true global optimum cannot be reached by gradient-based optimization. In their example, the true optimum can only be reached by removing one of the members. However, the stress constraint on that specific member prevents reducing the cross-sectional area of that member to zero. Kirsch (1989Kirsch ( , 1990 showed that singular optima are located in degenerate subspaces of the feasible domain, which are of a lower dimension than the design space. In density-based topology optimization (Bendsøe 1989), the presence of singular optima prevents the optimizer to reduce densities to zero; in general, large areas of intermediate densities will be present in the final design (Duysinx and Bendsøe 1998). We refer to Rozvany and Birker (1994) and Rozvany (2001) for extensive studies on singular optima and their fundamental characteristics.
A variety of solutions have been proposed to solve these fundamental difficulties. Currently, the conventional approach is to apply (i) relaxation to make singular optima accessible, and (ii) aggregation techniques to reduce the computational costs. Relaxation perturbs the feasible domain of the original optimization problem by replacing the original constraints by smooth approximations, which are always satisfied when material becomes void. Relaxation techniques that have been used are -relaxation (Cheng and Guo 1997), the qp-approach (Bruggi 2008), and defining a 'relaxed stress' (Le et al. 2010;Luo and Kang 2012).
Aggregation techniques reduce the number of constraints by lumping the local function values (stresses or constraint function values) into a single aggregation function. This aggregation function approximates the maximum local function value. The accuracy of this approximation depends on an aggregation parameter, and becomes exact in the limit as the aggregation parameter tends to infinity. In numerical practice, a moderate value is chosen for this aggregation parameter, which is a trade-off between two conflicting requirements: (i) an accurate approximation, and (ii) a sufficiently smooth function to prevent numerical difficulties. Aggregation functions that have been used are the P -norm (Duysinx and Sigmund 1998) and KS-function (Yang and Chen 1996). This transformation greatly reduces the number of constraints, and associated computational costs. However, since the aggregation parameter value is a trade-off, the aggregation function may not always be a sufficiently accurate approximation.
Recent, a lot of research has been aimed at improving the accuracy of aggregation functions. For example, using block constraints (París et al. 2010) in which the domain is subdivided into physical regions/blocks. Constraint aggregation is then applied on every subregion as a compromise between considering every local stress constraint and a single aggregation function. Le et al. (2010) and Holmberg et al. (2013) investigated similar approaches in which the composition of every region is based on the order of the stress values. Le et al. (2010) also proposed an adaptive normalization approach in which the aggregation function is scaled such that the maximum local stress converges to the allowable stress. Finally, Luo et al. (2013) proposed an enhanced aggregation method in which they combine an active set strategy with constraint aggregation. Using these techniques, improved optimized designs were obtained in which the maximum stress is close to the allowable stress. Nevertheless, the optimal settings such as the optimal number of regions may be very problem dependent and difficult to determine a priori. Furthermore, the computational costs increases with the number of regions. These difficulties have motivated the present work, which aims to explore a different direction to address stress-constrained topology optimization, thus providing an alternative approach with characteristics that differ from existing techniques.. This paper proposes a new method for topology optimization with local failure constraints. The general concept is to penalize the presence of local failure in a mechanical body by damaging material where local failure occurs. Damaging material here means locally degrading the material properties depending on the amount of local failure. We degrade material in an additional, so-called, damaged model of the same mechanical body. Assuming that the overall performance of the structure (e.g., compliance) is a monotonic function of the local material properties, degraded material will never improve its overall performance. Consequently, the damaged model will always perform worse than (or at best equally as) the original undamaged model. Following this idea, we can prevent local constraint violation indirectly by imposing a single constraint that both models should have the same overall performance. This constraint then prevents local failure in material regions that contribute to the overall performance of the structure. Therefore, we can obtain a minimum mass design that satisfies local constraints by minimizing mass under this condition of equal overall performance.
The concept of damage in the proposed method serves merely as a mechanism to penalize local constraint violation. We do not focus at accurately model a physical damage process since we aim for an optimized design without damage. Therefore, the proposed method is closer related to stress-based topology optimization methods, which aim at preventing yield failure than recent contributions in topology optimization considering nonlinear continuum damage mechanics (Amir and Sigmund 2013;James and Waisman 2014).
We apply the proposed damage approach to topology optimization with local stress constraints, and use the compliance as a measure of the overall performance. However, expectations are that the general concept can be applied to a wider range of problems with other local constraints, such as local temperature constraints in thermal problems. We validated the method on a three-bar truss example, and tested it on several benchmark problems in density-based topology optimization. Although a full comparative study with other approaches is outside the scope of this paper, our results show that the method is a viable alternative for conventional methodologies using relaxation followed by aggregation. The approach has the similar advantage that only a single performance constraint is considered instead of local stress constraints individually.
The remainder of this paper is structured as follows. Section 2 describes the proposed damage approach conceptually. Section 3 discusses numerical implementation aspects. In Section 4 we validate the method on a threebar truss optimization problem (Kirsch 1990). In Section 5 we discuss its implementation in density-based topology optimization, which includes the sensitivity analysis and the associated computational costs. In Section 6 the damage approach is applied on different numerical examples. Finally, conclusions are drawn in Section 7.

Damage approach in stress-constrained topology optimization
In this section, we discuss the damage approach applied to stress-constrained topology optimization. First, we discuss the concept of degrading material in which the stress exceeds the allowable stress. Next, we discuss how to formulate the optimization problem.

Damaged model
Consider a mechanical body of an isotropic elastic material that occupies the (design) domain ∈ R d (d = 2, or 3) with a boundary that consists of two disjoint parts: = D ∪ N . We define a traction force t on N , and a prescribed displacement on D . For simplicity, we assume the absence of body forces. Finally, we assume that material failure occurs once an equivalent stress criterion (e.g., Von Mises stress) exceeds the allowable stress: |σ (x)| > σ lim .
In the damage approach two different models describe the same mechanical body: the original model and the damaged model. The original model is shown in Fig. 1a in which E(x) denotes the strictly positive Young's modulus. In density-based topology optimization, this is the effective Young's modulus that may vary locally since it depends on an underlying density field: ρ(x) ∈ (0, 1]. We can then find a unique displacement field u that satisfies the boundary value problem associated with the original model. From this displacement field, we derive an equivalent stress criterion σ (x). The red region then denotes the subdomain where the stress exceeds the allowable stress: σ ⊆ . The next step is to degrade material in the regions where the stress exceeds the allowable stress. For this purpose we introduce the damaged model in Fig. 1b. All quantities associated with this damaged model have a tilde. In the damaged model, we define the Young's modulus such that it satisfies the following condition: Here,Ẽ is strictly positive, and smaller or equal than the Young's modulus in the original model E. For the damaged model we can now also set up boundary value problem and find a unique displacement field u. Notice that the original model remains undamaged; i.e., stress violation in the original model only affects the Young's moduli in the damaged model.
Suppose that the overall performance of the structure can be measured by a scalar function that depends monotonically on the local material properties. In that case, following (1), the damaged model will never perform better than the original model. In this purely mechanical problem, we use the compliance as measure of the overall performance, since it depends monotonically on the Young's moduli. Consequently, the damaged model will be always more (or at best equally) compliant: where C andC denote the compliance of the original and damaged model, respectively. Next we use this inequality to define the optimization problem.

Optimization problem
Our aim is to find the lightest design without violating any local stress constraints. We have seen from (2) and (1), that material where the stress exceeds the allowable stress, leads to a more (or at best equally) compliant damaged model. Consequently, we can enforce local stress constraints indirectly by a single equality constraint, which states that both models should have the same compliance. The optimization problem is defined as Here, V is the volume objective and s are the design variables in the design space S; e.g., densities in topology optimization, or cross-sectional areas in truss optimization. We assume here that the equilibrium equations are satisfied for every state of the design variables in the design space; therefore, (3) is written in its nested form without including the equilibrium equations of both models as equality constraints. Finally, h denotes the equality constraint, which is satisfied as long as the local stress constraints are satisfied in material regions that contribute to the overall compliance.

Implementation
This section briefly discusses the implementation aspects to solve the problem stated in (3) using standard gradientbased optimization.

Material degradation
We implement material degradation in (1) by establishing a relationship between the Young's modulus of the damaged model and the original undamaged model as Here, E min denotes a small positive number that acts as a lower bound on the Young's modulus to avoid singularity of the global stiffness matrix. Furthermore, β is the damage function introduced to degrade material as a function of the ratio of a scalar stress criterion and allowable stress: |σ |/σ lim . For simplicity, but without loss of generality, we assume that degradation is based on a single stress value per element; e.g., the axial stress in a truss element, or an equivalent stress criterion evaluated at the centroid in continuum finite elements.
The damage function β should be chosen such that (4) satisfies (1). In order to solve the problem using gradient-based optimization additional criteria are necessary. The damage function should be (i) at least first order differentiable, (ii) bounded asymptotically from below by zero to be consistent with physics, and (iii) monotonically decreasing once the stress exceeds its allowable limit. Many functions satisfy these criteria. Figure 2 shows the damage function used in this work, which is Here, α > 0 is the degradation parameter which controls the steepness of the damage function; i.e., the amount of damage relative to the stress level.
We emphasize that we do not intent to accurately model a physical damage process since we aim for a design without damage. In the present context, the damage function should be regarded rather as a penalty function to drive the solution towards a design without stress constraint violation.

Modified optimization problem
To solve the optimization problem in (3) numerically, we consider a slightly modified optimization problem. In gen- eral, topology optimization problems can be solved efficiently using sequential convex programming algorithms, such as MMA (Svanberg 1987) and CONLIN (Fleury 1989). The standard forms of these algorithms do not support nonlinear equality constraints directly. A solution is to replace the equality constraint in (3) by a pair of inequalities: h ≤ 0 and h ≥ 0. Since, the second inequality is true by definition (see (2)), the following equivalent problem can be formulated: Here, we have also introduced a small positive parameter δ to relax this inequality. By relaxing the constraint, the constraint is made less strict and inactive when there is no local stress constraint violation. Without relaxation, the constraint g(s) would be always active or violated.
In Section 4, we will show that this relaxation is necessary to make the optimum accessible for gradient-based optimization.

Validation on three-bar truss example
Before discussing the implementation of the method in density-based topology optimization, we consider the threebar truss problem shown in Fig. 3 (Kirsch 1990). We assume linear elasticity and small displacements. The general objective is to find the lightest structure in which the stress in none of the members exceeds the allowable stress. This example is a well-known benchmark in stress-constrained optimization since the true optimal solution is a so-called singular optimum, which cannot be accessed by gradientbased optimization. First, we consider the design space with the stress constraints as a reference, and show the effect of constraint relaxation on the feasible domain. Next, we investigate the feasible domain for the damage approach. Details on the calculation of stress and compliance are given in Appendix A.

Fig. 3
Three-bar truss introduced by Kirsch (1990) to demonstrate the existence of singular optima

Stress-constrained topology optimization
The objective is to minimize the mass m subject to the stress constraints g j : Here, d , denotes the set of all structural members in the discretized design domain. For every member, ρ e denotes the material density, A e the cross-sectional area and L e the length. The design variable vector A = (A 1 , A 2 ) contains the cross-sectional areas. Here, the crosssectional area of the third member and second member are the same: A 3 = A 2 . Finally, stress constraints g j are imposed only on the subset K ⊆ d of members with a strictly positive cross-sectional area. Since this subset depends on the current design A, the stressconstrained problem is known as an optimization problem with 'design-dependent constraints' (Rozvany 2001), also called 'vanishing constraints' (Achtziger and Kanzow 2008). In Fig. 3 the values are listed for the three-bar truss problem along with the Young's modulus E i for each member. Figure 4a shows the design space for the three-bar truss. The grey lines are the contour lines of the mass objective. The red, blue and green line are associated with the stress constraints in (7). The dashed lines indicate which side corresponds to constraint violation.

A singular optimum
Starting from an arbitrary point in the 'main' body of the feasible domain (the region above the first constraint g 1 ) a gradient-based optimizer typically converges to point D at A D = (0, 3/2), for which the mass is m D = 3/2. However, g 1 does not apply at this point since the cross-section A 1 becomes zero and the first member vanishes. In fact, the true optimum is found at point B located at A B = (0, √ 2/2), for which m B = √ 2/2. Point B is an example of what is known as a 'singular optimum' (Kirsch 1989(Kirsch , 1990). Therefore, the feasible domain for this problem consists of the region above the red line and the line segment B − D, as shown in Fig. 4b.
The difficulty of singular optima is that they are located in a degenerate subdomain of the feasible domain, such as line segment B − D in this example. These degenerate subdomains are of a lower dimension than the main body of the feasible domain, which makes them inaccessible using standard gradient-based optimization.

Stress constraint relaxation
The common way to handle singular optima is by constraint relaxation, e.g., -relaxation (Cheng and Guo 1997) and qp-relaxation (Bruggi 2008). The original stress constraints are replaced by approximations that are always satisfied for zero cross-section. The design space for such a relaxed problem does not contain any degenerate subdomains. Next, we will demonstrate the effect of −relaxation on the design space and the accessibility of the true optimum B in Fig. 4a.
For problems with vanishing constraints, one can reformulate the design-dependent set of stress constraints in (7) into the equivalent design-independent set of constraints (Cheng and Jiang 1992;Achtziger and Kanzow 2008): Here, the subset K has been eliminated and the new constraints g i are imposed on the entire set of structural members d . Every constraint g i is automatically  Fig. 4b. Therefore, the true optimum is still inaccessible to gradient-based optimization. However, now we can relax the constraints in (8) by introducing a positive relaxation parameter , which yields the −relaxed optimization problem: Notice that for every > 0, the constraint g ε j is always satisfied for a sufficiently small A j . The effect of relaxation is 'widening' the degenerate subdomains to the dimension of the design space. Figure 5 shows this effect for = 0.01 on the original constraints. The global optimum of the relaxed problem is close to the true optimum B and is accessible to gradient-based optimization.
This example demonstrates that constraint relaxation makes singular optima accessible. Cheng and Guo (1997) have demonstrated that the global optimum of the relaxed problem converges to the global optimum of the original problem (7) as approaches zero. However, Stolpe and Svanberg (2001) have shown that the trajectory of the global solution might be discontinuous. This trajectory is defined as the path of the global optimum to the relaxation parameter. As a consequence, finding the global optimum of the relaxed problem does not guarantee finding the global optimum of the original problem by following the path of this optimum as is decreased to zero.

Damage approach
Next, we apply the damage approach on this three-bar truss problem. Instead of three stress constraints, we only impose a single performance constraint: After substitution of the structural parameters listed in Fig. 3, the compliance of the original model and the damaged model are given by and respectively (see Appendix A). Here, the Young's moduliẼ i are damaged according to the damage function β defined in (5).
In the damage approach we can vary two parameters: the degradation parameter α > 0 and the relaxation parameter δ ≥ 0. Next, we investigate the effect of α and δ on the feasible domain and associated optima.

Design space of the unrelaxed problem: δ = 0
First, we consider the optimization problem in (10) without relaxation; i.e., δ = 0 for any α > 0. Figure 6 shows the design space for the damage approach. For δ = 0, the feasible domain and corresponding optima are independent of the value of α > 0 since no damage is allowed. The cyan color represents the region where the constraintg is active. Since there is no part of the design space where Fig. 6 Design space for three-bar truss problem for δ = 0 for any α > 0 the constraint is inactive, the cyan colored region represents the entire feasible domain. The former stress constraints are shown as dashed lines to indicate that they are not imposed directly in the problem formulation.
We observe that the feasible domain in the damage approach coincides with the feasible domain of the original optimization problem with local stress constraints (cf. Fig. 4b). However, here we only considered a single constraint. The feasible domain includes the line-segment B − D since stress violation there occurs in a 'vanished' member, and therefore, does not violate the constraint in (10). Since both feasible domains coincide, also in the damage approach, the optimum in general is a singular optimum (Kirsch 1989) inaccessible to standard gradient-based optimizers. Next, we demonstrate that relaxing the constraint by choosing δ > 0 makes singular optima accessible.

Constraint relaxation:
δ > 0 Figure 7a shows the effect of the relaxation parameter on the original constraint. The cyan colored lines represent the relaxed constraints for different values of δ > 0, and a fixed value of α = 1. The gray region represents the original feasible domain. We observe that relaxing the constraint widens subdomain B − D, and makes it accessible to gradient-based optimization. In contrast to the unrelaxed problem where the constraint surface (g = 0) occupies the entire feasible domain, the constraint surface is now a line (cf. Fig. 6), and all points above this line are feasible points where the constraint is inactive. As δ → 0, the perturbed feasible domain converges to the original feasible domain, which is true for any α > 0.
For the relaxed problem the found solution will contain a certain amount of damage since at the active constraint: C = C(1 + δ). This can observed in Fig. 7a as the relaxed constraints lie outside the original feasible domain. Consequently, optima for the relaxed problem will always violate a local stress constraint in at least one of the non-zero members. The amount of stress constraint violation depends on the degradation parameter α. For increasing values of α, material is damaged more rapidly once the stress constraint is violated. As a consequence, less stress violation is allowed. Figure 7b shows this effect for δ = 1 and increasing values of the degradation parameter α. We observe that the perturbed constraint becomes more conservative; i.e., less stress violation will be present in the optimized design.

A proper choice for α and δ
Our study on the effect of α and δ on the feasible domain and associated optima, demonstrated two effects: (i) the fea- sible domain of the relaxed optimization problem converges to the original feasible domain as δ → 0 for any fixed α > 0, and (ii) the relaxed constraint becomes more conservative as α increases for a fixed δ, which results in less stress constraint violation in the optimized design.
In this work, we exploit the first effect, and choose a small value of α. The reason is that for large values of α numerical instabilities may arise caused by large gradients along the constraint surface.

Density-based topology optimization
In this section, we present the method in the context of density-based topology optimization. First, we discuss the microscopic stress definition we have used. Finally, we discuss the sensitivity analysis and associated computational costs.

Density-based topology optimization
First, we establish a relationship between the stiffness of the discretized versions of the original model and the damaged model. For the original model in Fig. 1a, we adopt the modified SIMP interpolation scheme (Sigmund 2007). The design domain is discretized into finite elements and a density variable is assigned to each element, which can continuously vary between zero and one, representing 'void' and 'solid' material, respectively. The effective Young's modulus for each element in the original model is defined as Here, E 0 denotes the Young's modulus associated with solid material (ρ = 1), and E min is a small positive number that acts as a lower bound on the Young's modulus to avoid singularity of the global stiffness matrix. Finally, p > 1 is a penalization exponent introduced to penalize intermediate densities and promote a zero-one solution. The linear structural problem is then defined as Here, K denotes the global stiffness matrix, K e the element stiffness matrix associated with a Young's modulus of unity, and f the nodal load vector. One can now solve this problem for nodal displacements and calculate the stress.

Stress definition
A difficulty in density-based topology optimization is that the stress is non-uniquely defined for intermediate densities. Assuming that the density design variable in SIMP represents the effective stiffness of a porous microstructure (Bendsøe and Sigmund 2003), one can distinguish the stress at a macroscopic-and microscopic level (Duysinx and Bendsøe 1998). The macroscopic stress is based on the homogenized material properties of the microstructure: Here, ε e denotes the strain vector and C(E * e ) denotes the elasticity matrix based on the effective (homogenized) Young's modulus. For simplicity, we assume that the effective Young's modulus is defined following the traditional SIMP model: E * = ρ p E 0 . It turns out that the macroscopic stress is not suitable for stress-constrained topology optimization, since it cannot be used to predict failure at the microscopic level for intermediate densities (Duysinx and Bendsøe 1998). Furthermore, Le et al. (2010) demonstrated that using the macroscopic stress leads to an all-void design. Duysinx and Bendsøe (1998) have proposed a stress model that circumvents these problems that arise when using the macroscopic stress. Their stress model mimics the behavior of the microscopic stress (or "local" stress) in a layered composite, which is the stress experienced at the microscopic level. They have shown that in accordance with the stress behavior in such material, the microscopic stress in density-based topology optimization should be: (i) inversely proportional to the density and (ii) attain a finite value at zero density. This last condition follows from studying the asymptotic behavior of the microscopic stress in porous material as the density goes to zero. A definition consistent with condition (i) is Here, the value of the exponent q is chosen to satisfy the second condition (ii), which is only possible for q = p. Hence, the microscopic stress is defined as This definition of the microscopic stress is physically consistent for non-zero densities. However, since the microscopic stress remains finite at zero density, the feasible region contains degenerate subdomains, which causes singularity problems as discussed in truss topology optimization in Section 4. In the conventional approach in which stress constraints are considered directly, one typically overcomes these problems by relaxing the individual constraints; using for example, qp-relaxation (Bruggi 2008) or -relaxation (Duysinx and Bendsøe 1998). Another common approach, which has the same effect, is to consider a relaxed stress measure in which one enforces zero stress at zero density. For example, following (Le et al. 2010), the relaxed stress is based on the stress in (16) with the condition q < p. One can consider the difference between the exponents, qp = p − q as a measure of the amount of stress relaxation. For qp = 0, the relaxed stress becomes the microscopic stress. Unfortunately, the relaxed stress lacks a physical interpretation for intermediate densities.
In the damage approach, we use the physically consistent microscopic stress in (17) directly, and we only relax the constraint that states that the compliance of the damaged model should be less or equal to that of the original model. For clarity, we only plot this microscopic stress in 'material elements', defined as ρ ≥ 1/2. The reason is that the microscopic stress model, although physically consistent for intermediate densities, is non-zero in the voids. This difficulty is equivalent to the non-zero 'limiting stress' in truss optimization (Cheng and Jiang 1992); i.e., in truss optimization the stress converges to a finite value at zero cross-sectional area, which correspond to a member with a infinitesimal cross-sectional area. Consequently, the stress values are neglected in members that are (almost) vanished.

Damage model
The next step is to degrade material in the regions where the stress exceeds the allowable stress. We use the same discretization for both models. The Young's modulus for an element in the damaged model is defined as Here, the damage function β e (σ e ) is defined as in (5), and σ e is an equivalent stress criterion per element, based on the microscopic stress in (17).
Following (18) we construct the global stiffness matrix of the damaged model as and define the structural problem as whereũ denotes the vector of nodal displacements of the damaged model.

Sensitivity analysis
In the discrete density-based setting, the topology optimization problem is now defined as where V denotes the relative volume, V 0 the total volume of the design domain, and v e the volume of a finite element. The compliances of the original model and the damaged model are defined as C = f T u, andC = f Tũ , respectively. The sensitivity of the constraint with respect to a density design variable ρ e , is given by The second term in this equation contains the total derivative of the compliance of the original model. The compliance problem is self-adjoint (Bendsøe and Sigmund 2003), and the derivative is defined as dC dρ e = −pρ where u e is the vector with the nodal displacements of element e in the original model. Next, we derive the total derivative of the compliance of the damaged modelC in first term of (22). Here, we must take into account the relation between the material properties of the damaged model and the stresses in the original model. We consider an equivalent stress criterion based on the microscopic stress σ (u(ρ)), which only depends implicitly on the densities. We switch to index notation and summation convention. The indices run over the number of degrees of freedom. We calculate the sensitivities using the adjoint method. The first step is to add the equilibrium equations multiplied by an adjoint vector to the compliance: Next, we take the derivative with respect to a density design variable ρ e , and collect the terms containing the displacement sensitivities, which gives Next, the adjoints are chosen to eliminate the displacement sensitivities; i.e., the last two terms in (25) should become zero, which gives In the first term we made use of the fact that the first adjoint problem in (25) is self-adjoint: μ = −ũ. The second adjoint λ can be obtained by solving the adjoint problem: Here, symmetry of the global stiffness matrix is used, and z is the pseudo-load vector in which every component z k can be constructed as a summation over the elemental contributions: In this paper, we consider the Von Mises stress, based on the microscopic stress tensor in (17) evaluated at the centroid of each finite element. For a derivation of ∂σ e /∂u k we refer to Duysinx and Bendsøe (1998). In summary, the total derivative of the constraint function in (22) only requires the solution of the adjoint problem in (27). The two other adjoint problems are selfadjoint. Thus, the total computational costs are dominated by solving three systems of equations of the same size: the equilibrium equations of the original model and damaged model, and a single adjoint problem. The computational costs are comparable to conventional methodologies applying constraint aggregation over two regions assuming that no information is re-used of the factorized matrices.

Results and discussion
This section discusses the results obtained by testing the damage approach on different design cases shown in Fig. 8. First, we investigated the effect of the parameter settings on the optimized designs considering the cantilever in Section 6.1. Afterwards, the performance of the method under mesh refinement was investigated considering the L-bracket (Duysinx and Bendsøe 1998) in Section 6.2. Subsequently, Section 6.3 discusses the results obtained for a multiple load case (Le et al. 2010). Finally, Section 6.4 briefly discusses similarities and differences of the proposed method with conventional methodologies using relaxation and aggregation. Table 1 lists the general settings that apply to every example unless stated otherwise.

Cantilever
We investigated the parameter dependency of the optimized designs considering the cantilever shown in Fig. 8a. The design domain was discretized into a fixed finite element mesh of 100 by 50 equally sized quadrilaterals. The optimization problem was to minimize volume under an allowable stress of σ lim = 0.5. Fixed regular mesh in which every element has the same dimension: 1 × 1 Thickness 1 Young's Modulus E 0 = 1 Young's Modulus voids E min = 10 −9 E 0 Poisson's ratio ν = 0.3 Equivalent stress criterion Von Mises stress based on the microscopic stress in (17) evaluated at the centroid Distributed loads All loads are distributed over a length of 5

The effect of varying the relaxation parameter δ
First, the problem was solved for a fixed value of the degradation parameter α = 5, and different values of the relaxation parameter δ. Figure 9 shows all optimized designs, and the associated Von Mises stress distribution in the material elements. We observe that initially as the relaxation parameter decreases better performing stress-based designs are obtained in which the stress distributions are more uniform, and the maximum stress exceeds the allowable stress less. The reason is that lower values of δ allow less overstressed material in the optimized design because of the constraint C ≤ (1 + δ)C. However, decreasing the relaxation parameter also tends to lead to an increased number of iterations, and eventually the optimization process failed to converge to a black and white design (see Fig. 9f). Therefore, ideally, the relaxation parameter should be small, but provide sufficient relaxation to facilitate convergence.

Scaling both parameters by the same scaling factor
We investigated the effect on the optimized designs when scaling both the degradation-and relaxation parameter by the same scaling factor. It turns out that scaling both parameters by the same scaling factor gives equivalent designs over a large range of scaling factor values. To demonstrate this invariance we used the optimized design for (α 0 , δ 0 ) = (5, 1/64) in Fig. 9d as a reference, and solved the optimization problem again scaling both parameters by the scaling factor values a = 10 −2 , 10 −3 , and 10 −4 . Figure 10 shows that for a = 10 −2 , and 10 −3 equivalent designs were obtained as the reference design in Fig.  9d. For a = 10 −4 the optimized design performs less well, which indicates that this invariance under scaling in computational practice is limited to a certain range of scaling factor values. Nevertheless, these results indicate that one mainly searches a for a suitable ratio between the degradation-and relaxation parameter, rather than searching for suitable values for both parameters individually. The next subsection demonstrates that this invariance under scaling can be exploited to accelerate convergence by adjusting the optimizer settings.

Accelerate convergence
Given a certain ratio between both parameters that results in well-performing stress-based designs, it is advisable to choose the value of α small, and choose δ accordingly. The reason is that for lower values of α the constraint function is smoother (i.e., smaller gradients along the constraint surface), and one can accelerate convergence more by adjusting the optimizer settings, whereas for larger values of α this leads more quickly to numerical instabilities.
In conclusion, we have found that the relaxation parameter should be chosen as small as possible with respect to the degradation parameter, but provide sufficient relaxation to facilitate convergence. Furthermore, we have found scaling the degradation-and relaxation parameter by the same constant gives equivalent designs over a large range Cantilever optimized with relaxed optimizer settings to accelerate convergence: asyincr = 1.8, asydecr = 0.9, and an external move limit of 0.15. of scaling factor values. Consequently, for a given problem one searches for a suitable ratio between both parameters rather than searching for the best values of both parameters individually. Finally, for a certain ratio convergence can be accelerated by choosing α small and δ accordingly.

L-bracket: mesh refinement
The effect of mesh refinement was investigated using the Lbracket (Duysinx and Bendsøe 1998) in Fig. 8b. The design domain was discretized into a regular mesh of equally sized quadrilaterals. We used the conservative optimizer settings in Table 1, and the reference parameter settings are: (α 0 , δ 0 ) = (5, 1/64). The number of elements is defined as N = 14/25n 2 , where n denotes the number of elements along the longest edge in both vertical and horizontal direction. Starting from a reference mesh with n 0 = 100 (i.e., N 0 = 6400 elements), mesh refinement was applied considering multiples of n 0 ; i.e., n = 2n 0 , 3n 0 , and 4n 0 , which in terms of elements is a refinement of N = 4N 0 , 9N 0 , and 16N 0 . The allowable stress was set to σ lim = 1. Figure 12 shows the optimized designs and the associated stress distributions. The reference design in Fig. 12a for N 0 = 6400 has a rounded shape near the reentrant corner, which prevents a stress peak to occur as can be seen from the uniform stress distribution. It is observed that under mesh refinement, while using the same values for α and δ, the optimized design perform less well. Although the stress is mostly uniform distributed with a value close to the allowable stress, the optimized designs have a less or no rounded shape near the reentrant corner. Consequently, the maximum stress level increases under mesh refinement.
Our hypothesis is that this increased stress violation is caused by the necessary constraint relaxation. Because of constraint relaxation, the compliance of the damaged model is restricted by the following inequality:C ≤ (1 + δ)C. In other words, constraint relaxation allows the compliance of the damaged model to exceed the compliance of the original model by a maximum of δ × C. According to this fraction of the compliance of the original model a certain amount of degraded material (i.e., stress violation) is allowed in the optimized design. The contribution of the local stiffness of a single element to the overall stiffness decreases with mesh refinement. Therefore, for a finer mesh more stress violation is allowed in the same amount of elements, or alternatively, the same amount of stress violation is allowed in more elements. Therefore, it is expected that by decreasing δ inversely proportional to the increase in number of elements, the maximum stress should remain more or less constant.
We tested the aforementioned hypothesis by again applying mesh refinement, but now decreasing δ inversely proportional to the increase in the number of elements. Figure 13 shows the results under the same mesh refinement. In this case, the maximum stress increased less when refining the original mesh to 4N 0 , and is indeed nearly constant when refining the mesh from 4N 0 to 9N 0 , and 16N 0 . We observe that decreasing δ goes at the expense of an increased number of iterations. Slower convergence is caused by the fact that for smaller values of δ, the original feasible domain containing singular optima is perturbed less (see Fig. 7 for the truss example). This makes it more difficult for the optimizer to access correct optima and converge to a black and white design as was shown in Section 6.1.1.
In conclusion, the best ratio between the degradation and relaxation parameter is found to be mesh-dependent. The reason for this mesh dependency is that degradation Fig. 13 Mesh refinement applied to the L-bracket result in Fig. 12a while simultaneously decreasing δ. The parameter settings for the reference design are (α 0 , δ 0 ) = (5, 1/64), and N 0 = 6400 of a single element has less effect on the overall stiffness under mesh refinement. Therefore, under the same parameter settings, a finer mesh leads to an increased local stress violation. It was demonstrated that decreasing the relaxation parameter along with mesh refinement helps to maintain control over the maximum stress level. However, decreasing the relaxation parameter with respect to the degradation parameter also leads to an increased number of iterations.

Multiple load case
The optimization problem for multiple load cases is stated as where M is the number of constraints, which corresponds to the number of load cases. We solved the multiple load case in Fig. 8c (Le et al. 2010) for an allowable stress is σ lim = 1. The parameter settings were set to (α, δ) = 10 −3 ×(5, 1/64), and to accelerate convergence, the default optimizer settings were changed to asyincr = 1.8, asydecr = 0.9. In addition, the external move limit was changed to 0.15, and the convergence criterion was relaxed to: ρ ∞ < 0.05. Furthermore, the initial density field was set to ρ = 0.3 such that the initial design is infeasible.
In all previous results the maximum stress exceeded the allowable stress. In the next example, we applied an adaptive normalization strategy similar to Le et al. (2010) to converge closer to the allowable stress. Instead of degrading material where the stress exceeds the allowable stress a scaled version of the allowable stress is considered: Here, c is a scaling factor which is updated every x iterations by an increment c ∈ [− c max , c max ], which is defined as The increment is bounded from above and below by c max to avoid too large oscillations caused by changing the optimization problem.
For the multiple load case, we consider one scale factor for each load case. The initial scale factors were set to c 0 = 0.75, with a move limit of c max = 0.01. Adaptive scaling was initiated after first time convergence, and then the scale factors were changed adaptively every x = 3 iterations. Finally, to reduce oscillations a scale factor c i is only updated when |σ i max − σ lim |/σ lim > 10 −2 for that particular load case. Figure 14 shows the optimized design, and associated stress distribution for every load case. The optimized design has a rounded shape in both reentrant corners that prevent stress peaks to occur. For both load cases a uniform stress distribution was obtained. The maximum stress for both loading conditions was within 0.5 % of the allowable stress, which demonstrates that adaptive scaling strategies can be applied to converge close to the allowable stress. First time convergence was obtained after 103 iterations, and the total number of iterations was 117.

Damage approach vs. conventional approach
Although a detailed comparison with conventional methodologies is beyond the scope of this paper, we can discuss some similarities and differences between the proposed method and previously proposed methodologies based on constraint relaxation followed by constraint aggregation (e.g., Duysinx and Sigmund 1998;Le et al. 2010).
A similarity is that both approaches only consider a limited number of performance constraints, which drastically reduces the computational costs as opposed to considering every local stress constraint separately. The computational costs of the proposed approach are comparable to conventional methodologies applying constraint aggregation over two regions assuming that no information is re-used of the factorized matrices. In the proposed approach, constraint violation affects the compliance of the damaged model. Therefore, the compliance of the damaged model has a similar aggregation effect as conventional aggregation functions. A difficulty that arises in both methods is a loss of control over the maximum stress level as the number of local constraints increases.
As in previously proposed methodologies, constraint relaxation is necessary to make singular optima accessible. A difference is that in the damage approach constraint relaxation is a strictly mathematical procedure applied as a last step to be able to use gradient-based optimization. In conventional methodologies, relaxation is generally applied on the local constraints (or stresses) before applying constraint aggregation. For example, it has become common practice to consider a relaxed stress (Le et al. 2010;Luo and Kang 2012). In that case, relaxation is no longer a strictly mathematical procedure, but (slightly) alters the physics of the problem.

Conclusions
This paper presents a new method for topology optimization with local stress constraints. This method penalizes local stress constraint violation by degrading material where the stress exceeds the allowable stress. Since degraded material affects the overall compliance the optimizer promotes a design with a minimal amount of material degradation caused by local stress constraint violation. In this method, a limited number of performance constraints controls indirectly a large number of local stress constraints. Its behavior is controlled by two parameters that have a predictable effect on convergence and accuracy.
We validated the damage approach on an elementary truss example. It was demonstrated that the feasible domain of the damage approach coincides with that of the original stress-constrained topology optimization problem. Subsequently, the effectiveness of the method was demonstrated on several test problems in densitybased topology optimization, and results were obtained corresponding well with previously published stress-based results.
Optimized designs proved invariant to simultaneous scaling of the relaxation-and degradation parameter. Consequently, for a given problem, an optimal ratio of both parameters is sought rather than optimal values for both parameters individually. We demonstrated that this fact can be used to accelerate convergence. However, the optimal ratio of both parameters was found to be mesh dependent. For a finer mesh, decreasing the relaxation parameter is necessary to maintain control over the local stress level, however, this also leads to an increased number of iterations. How to maintain a high convergence rate and accurate control of local stresses under mesh refinement is a topic of future research.
Finally, in this paper we have applied the method to topology optimization problems with stress constraints, but we expect that the approach can be applied to a wider range of problems. In general, the concept should be applicable to problems with local constraints imposed on the material domain in which the overall performance of the design depends monotonically on some local material property.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

A Kirsch' three-bar truss
Using static equilibrium and the compatibility conditions at the joint, the internal forces are derived as Here, we already made use of the fact that the length of all members are equal (L i = L). The stresses in the original model, with substitution of the structural parameters, are then defined as Depending on these stresses, we degrade the Young's moduli in the damaged model:Ẽ i . The next step is compute the compliances for both models. We first derive the downward displacement v at the joint: With substitution of the structural parameters, the compliance of the original model is then and the compliance of the damaged model is Here,ṽ is the vertical displacement of the damaged model. The Young's moduli are degraded according to the stress state in the original model:Ẽ i = β(σ i )E i , where, β ∈ [0, 1].