Structural Topology Optimization with Local Finite-Life Fatigue Constraints

: To improve the fatigue resistance of engineering structures, topology optimization has always been an effective design strategy. The direct calculation of large-scale local fatigue constraints remains a challenge due to high computational cost. In the past, the constraint aggregation techniques, such as the P-norm method, were often applied to aggregate local fatigue constraints into a global constraint, whereas the resultant optimal solution was not consistent with the original problem. In order to meet the local fatigue constraints accurately and reduce the number of constraints, the augmented Lagrangian scheme is employed to transform the original problem into the unconstrained problem. To evaluate the fatigue strength at every material point of structures under the proportional load with variable amplitude, we adopt the Sines fatigue criterion based on the Palmgren–Miner linear damage assumption. In addition, we solve the fatigue-constrained topology optimization problem on the unstructured polygonal meshes, which are not sensitive to numerical instabilities, such as checkerboard patterns, compared with lower-order triangular and bilateral meshes. We provide some numerical examples to validate the potential of the presented method to solve the fatigue-constrained topology optimization problem. Numerical results demonstrate that the optimized designs considering local fatigue constraints have a higher ratio of fatigue resistance to material consumption than those obtained through the traditional P-norm method. Therefore, the proposed approach retaining the local nature of fatigue constraints is more beneﬁcial for realizing the efﬁcient material utilization in structural topology.


Introduction
The design space of engineering structures is greatly enriched by the structural topology optimization method for metal additive manufacturing, which takes into account the forming needs of complex structures and high-performance components. Moreover, this method is widely used in aviation, aerospace, transportation, nuclear power, and other fields [1,2]. At present, most of the topology optimization research work focuses on the maximization of structural stiffness with volume constraints. However, engineering structures often serve in the working environment of alternating loads. Generally, the design strategy with optimal stiffness cannot fully satisfy the fatigue resistance requirements, which seriously affects the reliability of the structure during its entire life cycle. Therefore, considering fatigue performance constraints in the lightweight design of structures has important scientific and engineering significance for enriching the strength design theory of engineering structures, and improving service performance. To reduce stress concentration and structural fracture caused by high stress values, topology optimization problems with stress constraints are important [3][4][5][6][7][8][9][10]. The difficulty of topology optimization design with stress constraints is that a large number of local constraints increase the computational cost of sensitivity analysis. In order to reduce the number of stress constraints, aggregative methods, such as P-norm [3][4][5] and K-S function [6], are generally used to convert local stress constraints into a global constraint. In order to preserve the local nature of stress constraints, Silva et al. [7,8] and Emmrndoerfer et al. [9] used the augmented Lagrangian (AL) method to transform the multi-constrained problem into a series of unconstrained subproblems. Recently, Oliver et al. [10] combined the unstructured polygon mesh method [11] with the augmented Lagrangian method, as a stress aggregation-free technique, and proposed a topology optimization method, which is applied to engineering structures with complex boundaries. Other few studies on local structural optimization can be found in [12][13][14] based on the AL method. Different from the stress-constrained topology optimization problems, the fatigue failure has a strong correlation with the loading history, and the local features are more significant, which causes fatigue constraints topology optimization problems to be highly nonlinear; therefore, there are relatively few related research works. Some scholars presented the P-norm based method for fatigue-constrained topology optimization problem employing various fatigue failure criteria, such as Palmgren-Miner [15][16][17], modified Goodman [18,19], Sines [20] and Murakami [21], considering initial defects through a quasi-static analysis method [22,23], frequency domain [24][25][26] or time domain equivalent static load, and other dynamic analysis methods [27,28] to predict the structural fatigue damage. Although the P-norm method can greatly reduce the computational cost of fatigue constraints and their sensitivity analysis, it fails to precisely control fatigue constraints.
Therefore, this paper proposes a topology optimization framework for local fatigue constraint problems based on the AL approach. Under the variable-amplitude proportional load, the fatigue damage of the structure is evaluated by the Palmgren-Miner fatigue criterion. Taking the material consumption as the objective function and the fatigue performance as the constraint condition, a lightweight design method for structures based on local fatigue constraints is proposed. In addition, the unstructured polygonal mesh method is integrated into the topology optimization model to realize the fatigue resistance topology optimization design of structures with complex geometric boundaries. Compared with the optimization results of the existing P-norm method, the effectiveness of the proposed aggregation-free approach for the fatigue constraints topology optimization problem is verified.

Optimization Model
In the framework of topology optimization proposed in this paper, with the goal of strong structural quantification, the aggregation-free local fatigue constraint problem (Q 1 ) and the P-norm aggregated global fatigue constraint problem (Q 2 ) are solved. A comparative analysis of the optimized solutions of these two problems verifies the effectiveness of the proposed method. where is the element area vector, 1 = [1] N is the unit column vector, N is the number of elements, and y(z), z is the element density variable and the design variable, respectively. m V (y(z)) is the volume interpolation function, g Dj (z, u) is the fatigue damage variable on the element j, the operator · P is the P-norm of the vector function number, and u is the node displacement vector.

Fatigue Failure Analysis
In order to estimate the high-cycle fatigue damage of the structure with variableamplitude proportional load, in this paper, the rainflow-counting method is used to determine the mean stress and stress amplitude with the multiaxial stress state, the Sines criterion is used to estimate the fatigue equivalent stress, and finally the Palmgren-Miner linear cumulative damage model is used to evaluate the fatigue failure of the structure. The high-cycle fatigue of engineering structures is generally manifested as low-stress failure. Thus, the linear quasi-static analysis method can be used in the research. The fatigue response of the structure under the variable-amplitude load is obtained by linearly superimposing the structural response with the reference load.
To develop the clear black-and-white design, we employ a volume interpolation function on the basis of the threshold projection approach [29], which is expressed as: m V (y(z)) = tanh(αη) + tanh(α(y(z) − η)) tanh(αη) + tanh(α(1 − η)) (2) where α is the threshold density and η controls the sharpness of the projection function around the threshold density. Additionally, the stiffness interpolation function is defined using SIMP function [10,11], which is formulated as: where m E uses an Ersatz number ε 1 to prohibit numerical instability when y(z) → 0 , p is the SIMP penalty factor.
In order to strengthen the well-posedness of the optimization problem, the design variable function is convolved with the linear 'hat' kernel function with the normalized mapping [30], so that the element density variable and the design variable have the following relationship: y = Pz (4) where P is the filter matrix, R is the filter radius, x − x k 2 is the Euclidean distance between the centroids of elements and k, and q is the filter index. In order to solve the structural topology optimization problem with complex design domain, the unstructured polygonal mesh method is used to discretize the spatial design domain, as shown in Figure 1. Talischi et al. [11] demonstrated that unstructured polygonal meshes are not sensitive to numerical instabilities, such as checkerboard patterns, compared with lower-order triangular and bilateral meshes. Particularly, non-uniform polygonal meshes can be generated by Lloyd's algorithm for discretization of design domain with complex boundaries. Oliver et al. [10], Senhora et al. [31], and Nguyen et al. [32] implemented local stress-constrained topology optimization on the basis of this polygonal discretization. Therefore, we adapt the polygonal composite finite element technique for topology optimization with fatigue constraints.    We consider the variable-amplitude load as shown in Figure 2, and use the quasistatic analysis method to solve the fatigue response of the structure with the reference load:  For the n-sided reference element, the Wachspress shape function at the node is defined as: where α i (ξ) is the interpolation function, ξ = ξ 1 ξ 2 is the interior point of the n-sided element, and A i (ξ), A i+1 (ξ) is the area of the corresponding shadow triangle (as shown in Figure 1a), which can be calculated according to Equation (9), where p i−1 = p 1,i−1 p 2,i−1 and p i = p 1,i p 2,i are the position vectors of the neighboring nodes of the n-sided element.
We consider the variable-amplitude load as shown in Figure 2, and use the quasi-static analysis method to solve the fatigue response of the structure with the reference load: where K(y(z)) is the global stiffness matrix and K e0 is the element stiffness matrix of the solid material.  In order to avoid the singular solution phenomenon of the optim modified qp relaxation technology is used to deal with the eleme elastic structure without prestress, the element reference stress und is expressed as: where D and B are the constitutive matrix and strain-displacem material, respectively, and u is the displacement matrix of the el In order to avoid the singular solution phenomenon of the optimization problem, the modified qp relaxation technology is used to deal with the element stress. For a linear Mathematics 2023, 11, 1220 5 of 17 elastic structure without prestress, the element reference stress under the reference load is expressed as:σ e = m E (y(z e ))DBu e (12) whereD and B are the constitutive matrix and strain-displacement matrix of the solid material, respectively, and u e is the displacement matrix of the element nodes. In order to predict the fatigue damage of engineering structures with variable-amplitude loads, the peak and valley values of the stress spectrum are often extracted by the rainflowcounting method (as shown in Figure 2), and the stress cycles of different cyclic characteristics are determined, thereby obtaining the stress amplitude scaling factor and average stress scaling factor. According to Equation (13), we estimate the stress state of any stress cycle.
where σ e a,i and σ e m,i are the stress amplitude and average stress of element e in the i-th cycle, respectively; c ai and c mi are the stress amplitude and average stress scaling factor of the i-th cycle, respectively. On the basis of the Sines finite-life fatigue criterion, the equivalent uniaxial stress of the finite-life in the plane stress state can be obtained from the alternating octahedral shear stress and the hydrostatic average stress [33].
where σ e i is the uniaxial equivalent stress of element e at the stress cycle i and β is the material parameter, which is recommended as 0.5. The relationship between uniaxial equivalent stress and fatigue life is described in the Basquin equation expressed as: where σ f is the fatigue strength coefficient, N e i is the fatigue life of the structure with stress cycle i, and b is the fatigue strength index. The Palmgren-Miner linear cumulative damage model can be used to assess the fatigue damage of each material point of the structure under the entire load spectrum. Correspondingly, the structure must meet the following fatigue constraints to avoid fatigue fracture.
where D e represents the cumulative damage of element e, c D represents the scaling parameter (generally greater than 1), and n i represents the number of cycles of load cycle i.

Fatigue Constraint
For the optimization problem Q 2 , the P-norm fatigue constraint can be expressed as: where P is the P-norm factor. A high P-norm factor is favorable for estimating the maximum value of the function, but it will also increase the nonlinearity of the optimization problem, which is difficult to solve due to convergence difficulties. In order to make the P-norm damage close to the maximum value of the actual damage, the adaptive constraint scaling technique [14,34] is generally used to correct the P-norm damage. Referring to the way Giraldo-Londono et al. [35] deal with local stress constraints, we construct polynomial local fatigue performance constraints on each finite element to avoid fatigue failure.
According to Equation (18), when the local fatigue constraint g Dj (z, u) 1, g j and Λ 3 j have the same order, which drives the optimizer to obtain a solution towards low damage. Similarly, when D j (z, u) → 1 , g j and Λ j have the same order, so that g j retains the characteristics of traditional linear constraint. To highlight these points, Figure 3 contrasts the conventional linear constraint with the cubic constraint, as illustrated in Equation (18). This figure demonstrates that the constraint violation domain is more significantly penalized by the cubic constraint than the traditional linear constraint. However, both the cubic and linear constraints function similarly to each other when Λ j → 0 .
solve due to convergence difficulties. In order to make the P-norm damage close to the maximum value of the actual damage, the adaptive constraint scaling technique [14,34] is generally used to correct the P-norm damage.
Referring to the way Giraldo-Londono et al. [35] deal with local stress constraints, we construct polynomial local fatigue performance constraints on each finite element to avoid fatigue failure. (18) According to Equation ( ,u j D z ,  j g and Λ j have the same order, so that  j g retains the characteristics of traditional linear constraint. To highlight these points, Figure 3 contrasts the conventional linear constraint with the cubic constraint, as illustrated in Equation (18). This figure demonstrates that the constraint violation domain is more significantly penalized by the cubic constraint than the traditional linear constraint. However, both the cubic and linear constraints function similarly to each other when Λ → 0 j . Figure 3. Comparison between traditional linear and cubic constraints.

Sensitivity Analysis
The gradient-based Globally Convergent Method of Moving Asymptotes algorithm (GCMMA) has good stability and robustness [36]. This work applies GCMMA to solve the fatigue-constrained topology optimization problem.
Since the optimization problem Q 1 is represented by density-related information through m V (y(z)) and m E (y(z)), it is convenient to formulate the sensitivity of objectivity function in the light of these fields. As such, we introduce the vector of element volume fractions V = m V (y(z)) and the vector of element stiffness parameters E = m E (y(z)) to obtain the sensitivity of the objective function by chain derivation rule as: where where J m E = diag(m E (y 1 ), · · · , m E (y N )) and J m V = diag m V (y 1 ), · · · , m V (y N ) . The sensitivity of polynomial fatigue constraints to design variables is expressed as: where We decompose the fatigue equivalent stress σ e i from Equation (14) into two parts: the equivalent stress amplitude σ a e i and the equivalent average stress σ m e i .
σ a e i = σ e ax,i − σ e yx,i 2 + σ 2 e ax,i + σ 2 e ay,i + 6τ 2 e a,i /2 (24a) By chain derivation rule, the partial derivative of D j (z, u) with respect to E is formulated as:

Augmented Lagrangian Method
As is well known, one of the main challenges to solving the fatigue-constrained topology optimization problem is the large number of local fatigue constrains imposed at each material point to avoid material failure. For the optimization problem Q 1 , to lessen the computational burden associated with local stress constraints, this work employs the augmented Lagrangian method. This method can obtain the optimal solution by solving a series of unconstrained optimization subproblems corresponding to the original constrained optimization problem. As a result, we search the minimum of AL function at the k-th step during optimization as follows [35,36]: where denotes the penalization term associated with local fatigue constraints. In the AL framework, the equality constraints corresponding to g j are given as: where The terms µ (k+1) and λ (k+1) j are the quadratic penalty factor and Lagrange multiplier, respectively. α > 1 is an update parameter, and µ max is a supremum used to prevent numerical instabilities. The GCMMA approach is used to solve the normalized unconstrained subproblem (30). The optimization terminates once the convergent criterion < Tol is satisfied. A schematic flowchart of the AL-based topology optimization framework is illustrated in Figure 4.
Read input Optimum topology Note that the penalization term ( ) ( ) , k P z u is normalized by the number of elements, namely the number of fatigue constraints N . This treatment attributes to the fact that a that the convergence is unfavorably influenced. Senhora et al. [31] demonstrated that this normalization can inhibit numerical instabilities due to excessively large N , and, subsequentially, Nguyen et al. [32] and Oliver et al. [10] efficiently solved the stress-constrained topology optimization problem based on this adjudgment. An analogous treatment to normalize the AL function is founded in da Silva et al. [8], where the objective function is amplified by multiplying the number of elements N . An alternative strategy for improving the robustness of the AL approach is proposed by da Silva et al. [7], where both the initial and maximal quadratic penalty factor are normalized regarding the number of elements. These works emphasize the significance of normalization treatment to enhance the robustness of the AL approach.
The sensitivity of the normalized AL function with respect to the element volume fractions and stiffness parameters can be achieved by the chain derivation rule and is represented as: where Note that the penalization term P (k) (z, u) is normalized by the number of elements, namely the number of fatigue constraints N. This treatment attributes to the fact that a considerably large N makes P (k) (z, u) dominant in the AL function L (k) (z, u), such that the convergence is unfavorably influenced. Senhora et al. [31] demonstrated that this normalization can inhibit numerical instabilities due to excessively large N, and, subsequentially, Nguyen et al. [32] and Oliver et al. [10] efficiently solved the stress-constrained topology optimization problem based on this adjudgment. An analogous treatment to normalize the AL function is founded in da Silva et al. [8], where the objective function is amplified by multiplying the number of elements N. An alternative strategy for improving the robustness of the AL approach is proposed by da Silva et al. [7], where both the initial and maximal quadratic penalty factor are normalized regarding the number of elements. These works emphasize the significance of normalization treatment to enhance the robustness of the AL approach.
The sensitivity of the normalized AL function with respect to the element volume fractions and stiffness parameters can be achieved by the chain derivation rule and is represented as: where From Equation (31), P (k) is independent from the volume fraction V , and thereby: To escape the costly calculation of ∂u/∂E , we adopt the adjoint method, whereupon the residual forms of the amplitude and mean equilibrium R a i , R m i are exploited.
where u a i , u m i denote the amplitude displacement and mean displacement under the reference load F r , respectively. F r,a i and F r,m i denote the amplitude and mean components of F r , respectively. Next, we introduce adjoint vectors corresponding to the amplitude and mean equilibrium, ξ a i and ξ m i , and rewrite Equation (36) as: where Since the computation of Lagrange multipliers are expensive in Equation (39), it is desirable to avoid directly solving the adjoint equation at every stress cycle i. Both the amplitude and mean displacements, u a i and u m i , are scaled by the displacement u under the reference load and can be obtained as demonstrated in Equation (13), i.e., Accordingly, we rewrite the adjoint term of Equation (39) as: To vanish the implicit terms containing ∂u/∂E from Equation (39), we assign the adjoint variable to satisfy the following equation: from which we attain the sensitivity of P (k) as: In terms of Equation (32), we have ∂h j (z, u)/∂u = 0 once g j (z, u) < −λ (k) j /µ (k) is met, and otherwise: ∂σ e a,i ∂u The terms ∂σ e a,i /∂u and ∂σ e m,i /∂u can be obtained from Equations (12) and (13), and expressed as:

Results and Discussion
The aggregation-free approach maintains the local nature of fatigue constraint, which ensures the local definition of fatigue constraint at each finite element, as stated in Q 1 . On the contrary, the constraint aggregation technique based the P-norm approach groups all the local fatigue constraints into a global constraint, which represents the whole design domain by roughly approximating the fatigue damage, as stated in Q 2 . As such, various strategies to fatigue constrains theoretically lead to the discrepancy in the optimum solution, which can be verified by the following numerical examples.

L-Bracket
As shown in Figure 5, the size parameter L of the L-shaped bracket is set to be 1 m. The material properties are as follows: elastic modulus E 0 = 70 GPa, Poisson's ratio µ = 0.3, fatigue strength coefficient σ f = 1350 MPa, fatigue strength exponent b = −0.086. A fixed constraint is applied to the upper end of the L-shaped bracket. This example includes a single load and multiple different loads acting at different locations of the structure, where the reference load F r is specified to be 20 kN. The fitter radius R is given as 0.065, the penalty parameter p as 4.5, and the projection parameter η ranges from one to ten, which is updated by increasing three at every five iterations.  (12) and (13), and expressed as:

Results and Discussion
The aggregation-free approach maintains the local nature of fatigue constraint, which ensures the local definition of fatigue constraint at each finite element, as stated in  1 .
On the contrary, the constraint aggregation technique based the P-norm approach groups all the local fatigue constraints into a global constraint, which represents the whole design domain by roughly approximating the fatigue damage, as stated in  2 . As such, various strategies to fatigue constrains theoretically lead to the discrepancy in the optimum solution, which can be verified by the following numerical examples.

L-Bracket
As shown in Figure 5, the size parameter L of the L-shaped bracket is set to be 1 m. The material properties are as follows: elastic modulus F is specified to be 20 kN. The fitter radius R is given as 0.065, the penalty parameter p as 4.5, and the projection parameter η ranges from one to ten, which is updated by increasing three at every five iterations. The design domain is meshed by 50,000 polygonal elements, and the transverse load is uniformly applied to the neighboring 8 elements to the loading application point to avoid stress concentration. k = 1000 random load given by  The design domain is meshed by 50,000 polygonal elements, and the transverse load is uniformly applied to the neighboring 8 elements to the loading application point to avoid stress concentration. k = 1000 random load given by F k = rand(F r , −F r ) is applied to the structure with the scaling parameter of load spectrum c D = 750. Figure 6 exhibits convergence histories of the objective functions for the optimization statements Q 1 and Q 2 with a P-norm factor of P = 8. Figure 7 shows the corresponding optimal topologies and normalized damage distribution of the structure subjected to a single load. These two optimum topologies differ from each other due to different strategies to fatigue constraints. Since the P-norm method increases the nonlinearity of fatigue constraint, the resultant topology is attained in more optimization iterations. In addition, the optimal topology obtained by aggregation-free approach is 21.3% lighter than that obtained by P-norm method. The former also manages to have a smaller portion of the design fully damaged in contrast to the latter. These results demonstrate that the optimum design considering local fatigue constraints presents a high ratio of the fatigue damage to the material amount, which verifies the effectiveness of the proposed approach.  2 with a P-norm factor of P = 8. Figure 7 shows the corresponding optimal topologies and normalized damage distribution of the structure subjected to a single load. These two optimum topologies differ from each other due to different strategies to fatigue constraints. Since the P-norm method increases the nonlinearity of fatigue constraint, the resultant topology is attained in more optimization iterations. In addition, the optimal topology obtained by aggregation-free approach is 21.3% lighter than that obtained by Pnorm method. The former also manages to have a smaller portion of the design fully damaged in contrast to the latter. These results demonstrate that the optimum design considering local fatigue constraints presents a high ratio of the fatigue damage to the material amount, which verifies the effectiveness of the proposed approach. To manifest the efficiency of this strategy for multiple loads, we perform the fatigueconstrained topology of the L-bracket subjected to multiple loads, as shown in Figure 5a. These results, as plotted in Figure 8, show that more beam-like members are located around the right part of the structure to counteract multiple different loads acting at the horizontal and vertical edges.   2 with a P-norm factor of P = 8. Figure 7 shows the corresponding optimal topologies and normalized damage distribution of the structure subjected to a single load. These two optimum topologies differ from each other due to different strategies to fatigue constraints. Since the P-norm method increases the nonlinearity of fatigue constraint, the resultant topology is attained in more optimization iterations. In addition, the optimal topology obtained by aggregation-free approach is 21.3% lighter than that obtained by Pnorm method. The former also manages to have a smaller portion of the design fully damaged in contrast to the latter. These results demonstrate that the optimum design considering local fatigue constraints presents a high ratio of the fatigue damage to the material amount, which verifies the effectiveness of the proposed approach. To manifest the efficiency of this strategy for multiple loads, we perform the fatigueconstrained topology of the L-bracket subjected to multiple loads, as shown in Figure 5a. These results, as plotted in Figure 8, show that more beam-like members are located around the right part of the structure to counteract multiple different loads acting at the horizontal and vertical edges.  To manifest the efficiency of this strategy for multiple loads, we perform the fatigueconstrained topology of the L-bracket subjected to multiple loads, as shown in Figure 5a. These results, as plotted in Figure 8, show that more beam-like members are located around the right part of the structure to counteract multiple different loads acting at the horizontal and vertical edges.  2 with a P-norm factor of P = 8. Figure 7 shows the corresponding optimal topologies and normalized damage distribution of the structure subjected to a single load. These two optimum topologies differ from each other due to different strategies to fatigue constraints. Since the P-norm method increases the nonlinearity of fatigue constraint, the resultant topology is attained in more optimization iterations. In addition, the optimal topology obtained by aggregation-free approach is 21.3% lighter than that obtained by Pnorm method. The former also manages to have a smaller portion of the design fully damaged in contrast to the latter. These results demonstrate that the optimum design considering local fatigue constraints presents a high ratio of the fatigue damage to the material amount, which verifies the effectiveness of the proposed approach. To manifest the efficiency of this strategy for multiple loads, we perform the fatigueconstrained topology of the L-bracket subjected to multiple loads, as shown in Figure 5a. These results, as plotted in Figure 8, show that more beam-like members are located around the right part of the structure to counteract multiple different loads acting at the horizontal and vertical edges.

Antenna Bracket
An antenna bracket with the geometrical parameters L = 4 m, H = 6 m and boundary conditions are depicted in Figure 9. The material properties are as follows: elastic modulus E 0 = 120 GPa, Poisson's ratio µ = 0.3, fatigue strength coefficient σ f = 1350 MPa, and fatigue strength exponent b = −0.086. The fitter radius R is defined as 0.22, the penalty parameter p as 3, and the projection parameter η ranges from two to ten, which is updated by increasing one at every five iterations. ary conditions are depicted in Figure 9. The material properties are as follows: elastic modulus = . The fitter radius R is defined as 0.22, the penalty parameter p as 3, and the projection parameter η ranges from two to ten, which is updated by increasing one at every five iterations. The design domain is discretized into N = 50,000 polygonal elements, and all the transverse loads are consistently applied to the eight elements. Here, we consider two loading conditions, as plotted in Figure 9. The random load is prescribed by . The P-norm factor is the same as in the first example. Figure 10 illustrates the optimum designs and normalized damage distribution for an antenna bracket subjected to a single load under the respective P-norm and local constraints. When compared with the P-norm method, the aggregation-free method generates a more efficient layout of the interior beam-like members, such that relatively less of full damage occurs in the entire design domain. Accordingly, the resulting optimum design has a more superior fatigue resistance. In addition, from the perspective of material consumption, the optimal design obtained by the aggregation-free method is 32.6% lighter than that obtained by the P-norm method. This further affirms that the presented approach is more efficient for the lightweight design of engineering structures. To examine the efficiency of the presented approach for multiple loading condition, we achieve the optimal topology and damage distribution of the structure under multiple loads. As depicted in Figure 11, the optimizer automatically strengthens the optimized design by increasing the beam-like members around the application points. The design domain is discretized into N = 50,000 polygonal elements, and all the transverse loads are consistently applied to the eight elements. Here, we consider two loading conditions, as plotted in Figure 9. The random load is prescribed by P k = rand(F r , −F r ) with k = 1000 and c D = 30, 000. The P-norm factor is the same as in the first example. Figure 10 illustrates the optimum designs and normalized damage distribution for an antenna bracket subjected to a single load under the respective P-norm and local constraints. When compared with the P-norm method, the aggregation-free method generates a more efficient layout of the interior beam-like members, such that relatively less of full damage occurs in the entire design domain. Accordingly, the resulting optimum design has a more superior fatigue resistance. In addition, from the perspective of material consumption, the optimal design obtained by the aggregation-free method is 32.6% lighter than that obtained by the P-norm method. This further affirms that the presented approach is more efficient for the lightweight design of engineering structures.
An antenna bracket with the geometrical parameters =4m L , =6m H and boundary conditions are depicted in Figure 9. The material properties are as follows: elastic modulus = . The fitter radius R is defined as 0.22, the penalty parameter p as 3, and the projection parameter η ranges from two to ten, which is updated by increasing one at every five iterations. The design domain is discretized into N = 50,000 polygonal elements, and all the transverse loads are consistently applied to the eight elements. Here, we consider two loading conditions, as plotted in Figure 9. The random load is prescribed by . The P-norm factor is the same as in the first example. Figure 10 illustrates the optimum designs and normalized damage distribution for an antenna bracket subjected to a single load under the respective P-norm and local constraints. When compared with the P-norm method, the aggregation-free method generates a more efficient layout of the interior beam-like members, such that relatively less of full damage occurs in the entire design domain. Accordingly, the resulting optimum design has a more superior fatigue resistance. In addition, from the perspective of material consumption, the optimal design obtained by the aggregation-free method is 32.6% lighter than that obtained by the P-norm method. This further affirms that the presented approach is more efficient for the lightweight design of engineering structures. To examine the efficiency of the presented approach for multiple loading condition, we achieve the optimal topology and damage distribution of the structure under multiple loads. As depicted in Figure 11, the optimizer automatically strengthens the optimized design by increasing the beam-like members around the application points. To examine the efficiency of the presented approach for multiple loading condition, we achieve the optimal topology and damage distribution of the structure under multiple loads. As depicted in Figure 11, the optimizer automatically strengthens the optimized design by increasing the beam-like members around the application points.
An antenna bracket with the geometrical parameters =4m L , =6m H and boundary conditions are depicted in Figure 9. The material properties are as follows: elastic modulus = . The fitter radius R is defined as 0.22, the penalty parameter p as 3, and the projection parameter η ranges from two to ten, which is updated by increasing one at every five iterations. The design domain is discretized into N = 50,000 polygonal elements, and all the transverse loads are consistently applied to the eight elements. Here, we consider two loading conditions, as plotted in Figure 9. The random load is prescribed by . The P-norm factor is the same as in the first example. Figure 10 illustrates the optimum designs and normalized damage distribution for an antenna bracket subjected to a single load under the respective P-norm and local constraints. When compared with the P-norm method, the aggregation-free method generates a more efficient layout of the interior beam-like members, such that relatively less of full damage occurs in the entire design domain. Accordingly, the resulting optimum design has a more superior fatigue resistance. In addition, from the perspective of material consumption, the optimal design obtained by the aggregation-free method is 32.6% lighter than that obtained by the P-norm method. This further affirms that the presented approach is more efficient for the lightweight design of engineering structures. To examine the efficiency of the presented approach for multiple loading condition, we achieve the optimal topology and damage distribution of the structure under multiple loads. As depicted in Figure 11, the optimizer automatically strengthens the optimized design by increasing the beam-like members around the application points.

Portal Frame
This benchmark example concerns a portal frame with L = 1.6 m, H = 0.8 m, and l = 0.15 m, and boundary conditions as plotted in Figure 12. The material properties are as follows: elastic modulus E 0 = 100 GPa, Poisson's ratio µ = 0.3, fatigue strength coefficient σ f = 1350 MPa, and fatigue strength exponent b = −0.086. This structure is fixed at the bottom and subjected to a transverse variable-amplitude load of F r = 70 kN at the center of the upper edge. The fitter radius R is defined as 0.046, the penalty parameter p as 4.5, and the projection parameter η ranges from one to ten, which is updated by increasing three at every five iterations.

Portal Frame
This benchmark example concerns a portal frame with =1. . Figure 13 compares the results obtained when exploiting either the P-norm method with P-norm factor P = 10 or the aggregation-free method as the fatigue constraint strategy. These two optimized designs have the stress-concentrated material efficiently deleted at the reentrant corner on the lower central part of the portal frame. In addition, it is obvious that these two designs are almost fully damaged. Nevertheless, the aggregation-free method removes the small slim substructures, which have a minor effect on the fatigue-induced damage. Particularly, the optimum design obtained for the aggregation-free method consumes 23.9% less of the material amount in contrast to that obtained for the P-norm method. As is known, the filter radius significantly affects the optimal results, and thus we examine the influence of a filter radius. Figure 14 illustrates the optimal topologies using various filter radii for the portal frame. As the filter radius increases, more elements are involved in the filter operation, which leads to different paths toward an optimized design. Relatively low filter radius produces some small features difficult to manufacture, as shown in Figure 14a,b, whereas high filter radius produces obscure topology with a relatively large value of the objective function. Consequently, it seems that the filter radius of 0.046 can achieve the desirable optimized design. The design domain is divided into N = 50,000 polygonal elements, and the transverse load is consistently distributed onto 12 elements. k = 1000 random loads are prescribed by P k = rand(F r , −F r ) with the load scaling parameter c D = 10, 000. Figure 13 compares the results obtained when exploiting either the P-norm method with P-norm factor P = 10 or the aggregation-free method as the fatigue constraint strategy. These two optimized designs have the stress-concentrated material efficiently deleted at the reentrant corner on the lower central part of the portal frame. In addition, it is obvious that these two designs are almost fully damaged. Nevertheless, the aggregation-free method removes the small slim substructures, which have a minor effect on the fatigue-induced damage. Particularly, the optimum design obtained for the aggregation-free method consumes 23.9% less of the material amount in contrast to that obtained for the P-norm method.

Portal Frame
This benchmark example concerns a portal frame with =1. . Figure 13 compares the results obtained when exploiting either the P-norm method with P-norm factor P = 10 or the aggregation-free method as the fatigue constraint strategy. These two optimized designs have the stress-concentrated material efficiently deleted at the reentrant corner on the lower central part of the portal frame. In addition, it is obvious that these two designs are almost fully damaged. Nevertheless, the aggregation-free method removes the small slim substructures, which have a minor effect on the fatigue-induced damage. Particularly, the optimum design obtained for the aggregation-free method consumes 23.9% less of the material amount in contrast to that obtained for the P-norm method. As is known, the filter radius significantly affects the optimal results, and thus we examine the influence of a filter radius. Figure 14 illustrates the optimal topologies using various filter radii for the portal frame. As the filter radius increases, more elements are involved in the filter operation, which leads to different paths toward an optimized design. Relatively low filter radius produces some small features difficult to manufacture, as shown in Figure 14a,b, whereas high filter radius produces obscure topology with a relatively large value of the objective function. Consequently, it seems that the filter radius of 0.046 can achieve the desirable optimized design. As is known, the filter radius significantly affects the optimal results, and thus we examine the influence of a filter radius. Figure 14 illustrates the optimal topologies using various filter radii for the portal frame. As the filter radius increases, more elements are involved in the filter operation, which leads to different paths toward an optimized design. Relatively low filter radius produces some small features difficult to manufacture, as shown in Figure 14a,b, whereas high filter radius produces obscure topology with a relatively large value of the objective function. Consequently, it seems that the filter radius of 0.046 can achieve the desirable optimized design.

Rectangular Plate with Round Hole
As depicted in Figure 15,

Rectangular Plate with Round Hole
As depicted in Figure 15, this example deals with a rectangular plate of L = 1.2 m and H = 0.4 m with circular holes of D = 0.24 m. This rectangular plate is pinned at the left edge and is subjected to a transverse variable-amplitude load at the center of the right free edge. The material properties are as follows: elastic modulus E 0 = 100 GPa, Poisson's ratio µ = 0.3, fatigue strength coefficient σ f = 1350 MPa, and fatigue strength index b = −0.086. The fitter radius R and the penalty parameter p are set to be 0.05 and 4.5, respectively. The projection parameter η varies from two to ten, which is updated by increasing one at every five iterations.

Rectangular Plate with Round Hole
As depicted in Figure 15 . The fitter radius R and the penalty parameter p are set to b respectively. The projection parameter η varies from two to ten, which increasing one at every five iterations. , where the reference load is r 15kN F = . Figure 16 shows the op and damage distribution for this rectangular plate with circular holes usin straint strategies. As is expected, the fatigue optimization attempts to elim damaged bar near the fixed boundary. These two designs resemble each oth and normalized damage distribution. The optimal topology obtained for th free method is 8.1% lighter than the counterpart obtained for the P-norm m mainly because the former configurates two additional declinate bars be circular holes, which strengthen the structure on the consumption of less ma

Curved Plate with Circular Holes
In this example, we consider a design of a curved plate with circula geometry and boundary are demonstrated in Figure 17  The design domain is divided into N = 30,000 elements, and the transverse load is imposed on the 6 elements around the center at the right vertical edge. k = 1000 random loads are applied by defining P k = rand(F r , −F r ) with the load scaling parameter c D = 10, 000, where the reference load is F r = 15 kN. Figure 16 shows the optimal topology and damage distribution for this rectangular plate with circular holes using various constraint strategies. As is expected, the fatigue optimization attempts to eliminate the undamaged bar near the fixed boundary. These two designs resemble each other in topology and normalized damage distribution. The optimal topology obtained for the aggregation-free method is 8.1% lighter than the counterpart obtained for the P-norm method. This is mainly because the former configurates two additional declinate bars between the two circular holes, which strengthen the structure on the consumption of less material amount.

Rectangular Plate with Round Hole
As depicted in Figure 15 . The fitter radius R and the penalty parameter p are set to be 0.05 and 4.5, respectively. The projection parameter η varies from two to ten, which is updated by increasing one at every five iterations. , where the reference load is r 15kN F = . Figure 16 shows the optimal topology and damage distribution for this rectangular plate with circular holes using various constraint strategies. As is expected, the fatigue optimization attempts to eliminate the undamaged bar near the fixed boundary. These two designs resemble each other in topology and normalized damage distribution. The optimal topology obtained for the aggregationfree method is 8.1% lighter than the counterpart obtained for the P-norm method. This is mainly because the former configurates two additional declinate bars between the two circular holes, which strengthen the structure on the consumption of less material amount.

Curved Plate with Circular Holes
In this example, we consider a design of a curved plate with circular holes, whose geometry and boundary are demonstrated in Figure 17

Curved Plate with Circular Holes
In this example, we consider a design of a curved plate with circular holes, whose geometry and boundary are demonstrated in Figure 17. The dimension parameters of the structure are as follows: L = 2.8 m, D 1 = 0.35 m, D 2 = 0.6 m, R 1 = 0.3 m, R 2 = 0.5 m. The material properties are as follows: elastic modulus E 0 = 100 GPa, Poisson's ratio µ = 0.3, fatigue strength coefficient σ f = 1350 MPa, fatigue strength index b = −0.086. This structure is clamped at the right circular hole and is subjected to a uniform transverse variable-amplitude load at the left circular hole. The fitter radius R and the penalty parameter p are specified as 0.17 and 3, respectively. The projection parameter η shifts from one to ten, which is updated by increasing three at every five iterations. This structure is clamped at the right circular hole and is subjected to a uniform transverse variable-amplitude load at the left circular hole. The fitter radius R and the penalty parameter p are specified as 0.17 and 3, respectively. The projection parameter η shifts from one to ten, which is updated by increasing three at every five iterations.  Figure 18 presents the optimal topology and damage distribution for the curved plate with circular holes using various constraint strategies. The optimized design obtained by the aggregation-free strategy is strengthened by adding several interior bars in external frames, where the optimized design by the P-norm method has high fatigue damage. Therefore, the former design is superior to the latter design in the fatigue resistance. Note that the aggregation-free strategy efficiently assigns the material towards impairing the fatigue damage, such that the resultant design consumes 21.8% less material than the counterpart obtained by the P-norm method.

Conclusions
A topology optimization framework is developed for local fatigue constraints employing the augmented Lagrangian (AL) approach. We introduce the AL approach into the optimization algorithm, whereupon the solution of the original constrained optimization problem is obtained by solving a series of unconstrained optimization subproblems. In contrast to the constraint aggregation technique, the aggregation-free strategy retains the local nature of fatigue constrains, such that the resulting optimum designs have slightly smaller portion of fatigue damage, while a considerable amount of material is saved, as demonstrated by numerical examples.
This optimization scheme based on the AL method is generic in Engineering design, where the stress-based fatigue criteria are employed. Note that the presented approach does not implement much more work than the constraint-aggregation-based fatigue optimization, and provides a powerful tool in conceptual design for engineering problems driven by fatigue.   The design domain is discretized into N = 30,000 polygonal elements. The structure is subjected to k = 1000 random loads defined as P k = rand(F r , −F r ), where the load scaling parameter is c D = 11, 000. The reference load is F r = 10 kN, which is uniformly applied onto the upper part of the left circular hole. Figure 18 presents the optimal topology and damage distribution for the curved plate with circular holes using various constraint strategies. The optimized design obtained by the aggregation-free strategy is strengthened by adding several interior bars in external frames, where the optimized design by the P-norm method has high fatigue damage. Therefore, the former design is superior to the latter design in the fatigue resistance. Note that the aggregation-free strategy efficiently assigns the material towards impairing the fatigue damage, such that the resultant design consumes 21.8% less material than the counterpart obtained by the P-norm method.
This structure is clamped at the right circular hole and is subjected to a uniform transverse variable-amplitude load at the left circular hole. The fitter radius R and the penalty parameter p are specified as 0.17 and 3, respectively. The projection parameter η shifts from one to ten, which is updated by increasing three at every five iterations. , which is uniformly applied onto the upper part of the left circular hole. Figure 18 presents the optimal topology and damage distribution for the curved plate with circular holes using various constraint strategies. The optimized design obtained by the aggregation-free strategy is strengthened by adding several interior bars in external frames, where the optimized design by the P-norm method has high fatigue damage. Therefore, the former design is superior to the latter design in the fatigue resistance. Note that the aggregation-free strategy efficiently assigns the material towards impairing the fatigue damage, such that the resultant design consumes 21.8% less material than the counterpart obtained by the P-norm method.

Conclusions
A topology optimization framework is developed for local fatigue constraints employing the augmented Lagrangian (AL) approach. We introduce the AL approach into the optimization algorithm, whereupon the solution of the original constrained optimization problem is obtained by solving a series of unconstrained optimization subproblems. In contrast to the constraint aggregation technique, the aggregation-free strategy retains the local nature of fatigue constrains, such that the resulting optimum designs have slightly smaller portion of fatigue damage, while a considerable amount of material is saved, as demonstrated by numerical examples.
This optimization scheme based on the AL method is generic in Engineering design, where the stress-based fatigue criteria are employed. Note that the presented approach does not implement much more work than the constraint-aggregation-based fatigue optimization, and provides a powerful tool in conceptual design for engineering problems driven by fatigue.

Conclusions
A topology optimization framework is developed for local fatigue constraints employing the augmented Lagrangian (AL) approach. We introduce the AL approach into the optimization algorithm, whereupon the solution of the original constrained optimization problem is obtained by solving a series of unconstrained optimization subproblems. In contrast to the constraint aggregation technique, the aggregation-free strategy retains the local nature of fatigue constrains, such that the resulting optimum designs have slightly smaller portion of fatigue damage, while a considerable amount of material is saved, as demonstrated by numerical examples.
This optimization scheme based on the AL method is generic in Engineering design, where the stress-based fatigue criteria are employed. Note that the presented approach does not implement much more work than the constraint-aggregation-based fatigue optimization, and provides a powerful tool in conceptual design for engineering problems driven by fatigue.
Funding: This research was funded by National Natural Science Foundation of China (52175502, 51505096) and Heilongjiang Natural Science Fund Joint Guidance Project (LH2020E064).

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to lab privacy.