Computer Modeling in Engineering & Sciences

DOI: 10.32604/cmes.2022.017630

ARTICLE

Topology Optimization with Aperiodic Load Fatigue Constraints Based on Bidirectional Evolutionary Structural Optimization

Yongxin Li1, Guoyun Zhou1, Tao Chang1,*, Liming Yang2 and Fenghe Wu1

1Heavy-Duty Intelligent Manufacturing Equipment Innovation Center of Hebei Province, School of Mechanical Engineering, Yanshan University, Qinhuangdao, China
2Space Star Technology Co., Ltd., Beijing, China
*Corresponding Author: Tao Chang. Email: changtao20210901@163.com
Received: 25 May 2021; Accepted: 03 August 2021

1  Introduction

For air vehicles, components or supports are exposed to aperiodic loads, and there are strict requirements for minimum service life. The weight of aircraft parts directly affects carrying travel, and lightweight research can greatly improve the power-to-weight capacity. Therefore, this paper studies the topology optimization problem with fatigue life as the constraint and the minimum volume/weight as the objective under aperiodic loading.

The field of topology optimization mainly focuses on minimizing structural compliance [1–3] to increase structural stiffness. At present, there are roughly three ways to solve topology optimization with fatigue constraints. First, fatigue life is directly used as an optimal parameter, which is usually realized by a heuristic algorithm because fatigue life is not differentiable from variables. For example, Haiba et al. [4] used the ESO method in 2005 to establish an element deletion criterion based on fatigue life to gradually remove high-life elements. In the second way, fatigue cumulative damage is set as the design parameter instead of fatigue life. Jacob et al. [5] used fatigue cumulative damage in 2016 to describe the fatigue response and P-norm aggregation function to agglomerate fatigue constraints on the variable density method. In 2018, Seung et al. [6] developed topology optimization of fatigue constraints under variable amplitude cyclic loads based on the equivalent static load method. In 2020, Chen et al. [7] penalized cumulative damage to achieve anti-fatigue topology optimization. However, cumulative damage is highly nonlinear with stress, which seriously affects optimization stability. In the third way, fatigue is equivalent to stress using the fatigue failure criterion. In 2014, Holmberg et al. [8] developed a variable density method taking critical fatigue stress and static stress as constraint conditions simultaneously. In 2015, Jeong et al. [9] proposed a topology optimization method considering dynamic fatigue and static failure criteria using average stress to represent fatigue based on the variable density method. In 2019, Nabaki et al. [10] used a modified Goodman failure criterion to study the maximum stiffness problem under volume and fatigue life constraints. In 2020, Zhao et al. [11] discussed topology optimization under a high cycle fatigue constraint based on the Crossland criterion, and the fatigue constraint was equivalent to the peak stress not exceeding the threshold.

However, the aperiodicity and randomness of the aperiodic load determine that its fatigue cannot be described by stress or strain directly and can only be counted by the rain-flow method and expressed by fatigue damage. Due to the nonlinearity of fatigue damage with life history, there are few studies concerning topology optimization under aperiodic load fatigue constraints. This research aims to address the high cycle fatigue problem under aperiodic loading. If the stress value can be connected with fatigue damage [8], the fatigue constraint can be transformed to a critical stress problem.

Stress constraint topology optimization has conundrums of locality. The stress of each element must be constrained to search for the maximum stress. The sensitivity of each constraint needs to be calculated, so it leads to excessive computation whether using a direct method or adjoint method. To solve this problem, aggregation functions, such as the P-norm method [12] and K-S function method [13], can be applied. In 2013, Luo et al. [14] proposed an enhanced K-S function based on the variable density method to solve topological optimization for minimizing volume under stress constraints. In 2019, Fan et al. [15] addressed minimizing compliance under volume and stress constraints based on the BESO method and P-norm equation. In 2019, Long et al. [16] studied stress-constrained topological optimization under a simple harmonic load using the P-norm function based on the variable density method. In 2020, Zhao et al. [17] improved BESO with a dynamic stress response under random loads by the P-norm function.

In this paper, aiming at the severe nonlinearity of cumulative damage under aperiodic loading, a transformational relation of fatigue damage with stress is established. First, fatigue damage of the substep is calculated by the stress time history, and the problem is transformed into a stress topology optimization problem based on BESO. In the classic BESO method, volume is generally the constraint, so volume evolution is adopted to realize element addition and deletion, and the objective change value is taken as the convergence criterion. However, volume is the objective here, and if volume is taken to measure evolution and convergence simultaneously, it will be difficult to converge due to missing stress. Therefore, a normalized parameter, which comprehensively considers volume and stress, is proposed and taken as the convergence parameter. In this way, the topology lightening problem with a high cycle fatigue life as a constraint is solved.

2  Fatigue Constraint Treatment

2.1 Fatigue Life Calculation under Aperiodic Load

To interpret the relationship between fatigue life, damage and stress, this section briefly introduces a fatigue analysis program for aperiodic loads. The general process is as follows: first, fatigue life is calculated by alternating stress and average stress extracted from the transient stress time history. To explain the multiaxial stress of each node, the signed von Mises stress is calculated. The reversal of symbolic von Mises stress and average stress are determined by rain flow counting. The Goodman correction method is used to calculate the amplitude of the effective stress. Then, the Basquin equation is applied to calculate the number of failure cycles for each effective stress amplitude. Finally, all the damage caused by stress reversal in the loading history is accumulated linearly by the Palmgren–Miner law [18].

In detail, the stress history is solved by a given dynamic system as

Mu¨(t)+Cu˙(t)+Ku(t)=F(t),(1)

where M, C, K, F(t), u(t) and t represent the uniform global mass matrix, global damping matrix, global stiffness matrix, time-varying force vector, displacement vector and time, respectively. The first and second derivatives of the displacement vector over time are expressed as u˙(t) and u¨(t). The initial values for displacement, velocity and acceleration are all zero. The implicit direct integration method is used here to solve the dynamic equation, which calculates the response history with step-by-step iterations. The displacement response of step n + 1 can be obtained and written as

K¯un+1=F¯n+1,(2)

where K¯ and F¯n + 1 are the effective stiffness matrix and effective vector of Step n + 1, respectively. Then, the stress vector of the i-th element of the n + 1 step σin+1 is calculated. Thus

σin+1=D0Buin+1(3)

where D0 and B represent the elastic matrix and strain matrix, respectively.

Equivalent parameters similar to von Mises stress are usually positive [19], so a negative sign of von Mises stress is determined by a positive or negative sign of the first invariant. The time history of alternating stress is obtained as

σsvm={σ1|σ1|⋅12(σ1−σ2)2+(σ3−σ2)2+(σ1−σ3)2if |σ1|>|σ3| σ3|σ3|⋅12(σ1−σ2)2+(σ3−σ2)2+(σ1−σ3)2if |σ3|>|σ1|,(4)

where σsvm, σ1, σ2 and σ3 are the von Mises stress and the first, second and third principal stresses, respectively. Stress amplitude σa and mean stress σm are determined by rain flow counting such that

σa=σmax−σmin2,(5)

σm=σmax+σmin2.(6)

To consider the influence of the average stress, the Goodman correction equation is used to calculate the amplitude of the effective stress

σa=σ−1(1−σmσb).(7)

where σ−1 is the constant life fatigue strength under symmetrical cycles and σb is the tensile strength.

The number of failure cycles (Ni) of each effective stress amplitude can be computed by the Basquin equation. Thus,

σa,i=σf(2Ni)b,(8)

where σf is the fatigue strength and b is the test constant. The failure times of the i-th cycle load Ni can be obtained as follows:

Ni=12(σa,iσf)1b.(9)

Cumulative fatigue damage can be calculated by the Palmgren–Miner law as

D=∑i=1lniNi.(10)

where l is load series and ni is the number of load cycles.

Therefore, fatigue damage is expressed by the stress amplitude as

D=∑i=1k2ni(σa,iσf)−1b.(11)

From Eq. (11), there is an exponential relationship between fatigue damage and stress, and it is nonlinear. Therefore, it is not conducive to setting fatigue damage as a constraint condition.

2.2 Fatigue to Stress Constraint

In this section, the process course of converting the fatigue lifetime requirement to a stress constraint is presented. First, the position of the danger point is obtained by transient dynamic analysis loading aperiodic cyclic force. In addition, the stress time history is obtained. Then, the rain flow counting method is used to extract the effective stress cycle. Fatigue damage is calculated with the Palmgren–Miner criterion. Finally, according to the lifetime requirements, the critical stress level is inversely calculated and taken as the constraint condition for topology optimization. The specific implementation steps are as follows:

(1)   With the requirements of fatigue design life, stress amplitude σa and cycle number n are derived linearly in the whole life cycle, and the maximum value of stress amplitude is found, which is recorded as S(0);

(2)   Assuming that k damage calculations are made, based on the Miner theory, the total damage can be obtained as D(0)=∑i=1k1Ni;

(3)   If D(0)<1, fatigue failure will not occur, so S(0) is increases and, vice versa, S(0) decreases;

(4)   S(0) takes a series of values S(1),S(2),…,S(m), the stress amplitudes change sequentially, and the corresponding damage D(1),D(2),…,D(m) is calculated until D(m)>1;

(5)   (D(i),S(i))(i=0,1,…,m) is fitted to a stress-damage curve, as shown Fig. 1. When D = 1, the corresponding S is predicted as the maximum stress of critical fatigue damage. Stress values are taken as constraint conditions, that is, the transformation from fatigue to stress constraint is realized as shown in Fig. 2.

Figure 1: Stress damage curve

Figure 2: Flow chart of transforming fatigue constraint to stress constraint

3  Fatigue Constrained Topology Optimization on BESO

Through the conversion above, fatigue-constrained topology optimization can be transformed into a stress-constrained problem. To reduce the quantity of stress constraints at every cycle, the P-norm aggregation function is used as

σPN=(∑j=1nσvm,jp)1p.(12)

where σvm,j is the Von Mises stress of the i-th element and p is the aggregation norm. Therefore, the mathematical formula of fatigue constrained topology optimization on the BESO method reads

find:xi=[x1,x2,⋯,xn]min:V=∑i=1nVixis.t.:M\textit{{\"U}}+CU+KU=F,σPN≤σvm∗xi=0 or 1(13)

where xi and Vi represent the i-th element variable and the i-th element volume respectively. σvm∗ is the critical stress after equivalence.

Classic BESO uses the change of objective function as the convergence criterion. The objective parameter of optimization model in this paper is volume. If the volume change is set as the convergence criterion, the volume increment between adjacent iterations will alter drastically, and the convergence value is hard to realize. To take into account both optimization objective and constraint parameters, this paper proposes a normalized performance index PI, which can comprehensively consider stress change and volume change as

PI=σPN0V0σPNiVi,(14)

where σPN0 and σPNi represent the initial and i-th stress aggregation values, respectively. V0 and Vi represent the volume in the initial and i-th iterations, respectively. Thus, the termination criterion of the optimization process reads

error=|∑i=1MPIk−i+1−∑i=1MPIk−M−i+1|∑i=1MPIk−i+1≤τ.(15)

In general, M refers to the range of iteration intervals that determines convergence, usually 3–10, which means that the performance change over the past 10 iterations is small enough to achieve convergence when it is less than τ. τ is the convergence error.

Except for sensitivity analysis, convergence criterion and fatigue constraint, the rest of the flow is based on BESO theory [20]. A flow chart of the optimization algorithm with fatigue constraints is shown in Fig. 3.

Figure 3: Flow chart of optimization algorithm with fatigue constraints

4  Numerical Examples

To show the effectiveness and correctness of the proposed method, we design two numerical examples: double-bar plates and L-shaped beams. Therefore, a double-bar plate is designed with no obvious stress concentration, and an L-shaped beam is designed with stress concentration. The fatigue life is designed to be 140 h to confine to elastic deformation. The stress aggregation norm in Eq. (12) is set at p = 6.

4.1 Double-Bar Plate

The initial design area and boundary conditions of the double-bar plate are shown in Fig. 4. The thickness is 1 mm. The structure is divided into 80 × 40 quadrilateral plane stress elements. A fixed constraint is applied on the lower side, and an aperiodic load is applied to the 11th node in the middle of the upper side. The time history of the aperiodic load is shown in Fig. 5. Young's modulus, Poisson's ratio and the material density are set as E = 71 GPa, v = 0.33 and ρ = 2793 Kg/m³. Material S-N curve equation expresses as

log⁡σa=8.951−0.103log⁡N.(16)

Figure 4: Initial design area of double-bar structure

Figure 5: Load time history

There is no stress concentration caused by the structure shape in the double-bar model. The iterative changes in the volume fraction and PI are shown in Fig. 6. After 100 iterations, the volume fraction and PI reach a platform, and the convergence criterion is met. When PI reaches approximately 3.36, it tends to converge, and the objective function is 0.23. Relatively, the maximum stress and P-norm aggregation stress show slight fluctuations during the last 20 iteration cycles, as shown in Fig. 7. The maximum stress value exhibits a smaller downward trend than the P-norm aggregation stress and an upward trend. This shows that the P-norm aggregation function helps local stress reach convergence. Fig. 8 shows the final optimal configuration and stress distribution result. Direction inclination leads to loading asymmetry between the left and right parts, resulting in the asymmetric configuration, as shown in Fig. 8a.

Figure 6: Variation curve of objective volume and PI

Figure 7: Variation curve of structural maximum stress and P-norm aggregation stress

Figure 8: Final iterative results of double-bar plate and stress nephogram (a) Final iterative configuration (b) Stress nephogram

4.2 L-Shaped Beam

The initial design area and boundary conditions of the L-shaped beam are shown in Fig. 9. The design area is meshed by four-node plane elements with a length of 1 mm. The upper side is fixed, and the aperiodic load is shown in Fig. 10. Force is applied on the top vertex of right boundary. To avoid stress concentration at the loading position, the load is evenly distributed to 8 nodes. The material parameters are exactly the same as those in the double-bar plate.

Figure 9: L-shaped beam structure

Figure 10: Load time curve

The variation curve of the volume fraction and PI is obtained, as shown in Fig. 11. With iteration progress, the volume and PI fluctuate violently after 30 steps. This is due to the high nonlinearity of stress. With optimization, in some areas with high stress concentrations, the high stress distribution is very sensitive to even minor topological changes. However, the amplitude of oscillation tends to decrease with increasing iteration number. When PI reaches approximately 1.73, it tends to converge, and the objective function is 0.69. Approximately 250 steps are consumed to reach convergence. Since the convergence criterion is a relative change between the last and first range of the iteration interval, iteration convergence occurs since a small floating value is reached. After repeated trial calculations, a better convergence effect is obtained when the range of iteration intervals M = 4 and convergence error τ = 0.0001.

Figure 11: Variation curve of objective function and performance index

Curves of the maximum stress and P-norm aggregation stress are shown in Fig. 12. The final topology and its stress cloud diagram are shown in Fig. 13. The maximum stress is 175.74 MPa, and the P-norm aggregation stress is 228.14.

Figure 12: Variation curve of structural maximum stress and P-norm aggregation stress

Figure 13: Final iterative result of L-shaped beam and stress nephogram (a) Final iterative shape (b) Stress nephogram

5  Conclusion

Fatigue life, strength and stiffness are three important factors that designers cannot ignore in engineering. Aiming at high cycle fatigue caused by aperiodic load, this paper realizes topology optimization method of aperiodic load based on BESO.

In order to solve the severe nonlinearity of fatigue damage to design variables, this paper uses fatigue analysis to calculate cumulative damage and then equivalent fatigue damage to critical stress by stress damage curve to transform fatigue constraint problem to stress constraint problem.

Then, P-norm is used for stress aggregation, and the dimensionless index PI is proposed to improve convergence speed. The effectiveness of proposed method in obtaining the lightest configuration which strictly satisfies life constraint is proved by two examples with and without stress concentration.

Funding Statement: This research is funded by Chinese National Natural Science Foundation (No. 51890881). Science and Technology Project of Hebei Education Department (Nos. ZD2020156, QN2018228).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.